Exponenta Banner

Уравнение Пуассона на единичном диске

Скрипт Open Live Script

В этом примере показано, как численно решить уравнение Пуассона, сравните числовое решение с точным решением и совершенствуйте mesh, пока решения не близки.

Уравнение Poisson на единичном диске с нулем граничное условие Дирихле может быть записано как -Δu=1 \in Ω, u=0 on δΩ, где Ω единичный диск. Точное решение

u(x,y)=1-x2-y24.

Для большинства УЧП не известно точное решение. Однако уравнение Пуассона на единичном диске имеет известное, точное решение, которое можно использовать, чтобы видеть, как ошибка уменьшается, когда вы совершенствовали mesh.

Описание задачи

Создайте модель PDE и включайте геометрию.

model = createpde();
geometryFromEdges(model,@circleg);

Постройте геометрию и отобразите метки ребра для использования в определении граничного условия.

figure 
pdegplot(model,'EdgeLabels','on'); 
axis equal

Figure contains an axes object. The axes object contains 5 objects of type line, text.

Задайте нуль граничные условия Дирихле на всех ребрах.

applyBoundaryCondition(model,'dirichlet', ...
                             'Edge',1:model.Geometry.NumEdges, ...
                             'u',0);

Задайте коэффициенты.

specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',1);

Решение и ошибка с крупной сеткой

Создайте mesh с целевым максимальным размером элемента 0.1.

hmax = 0.1;
generateMesh(model,'Hmax',hmax);
figure
pdemesh(model); 
axis equal

Figure contains an axes object. The axes object contains 2 objects of type line.

Решите УЧП и постройте решение.

results = solvepde(model);
u = results.NodalSolution;
pdeplot(model,'XYData',u)
title('Numerical Solution');
xlabel('x')
ylabel('y')

Figure contains an axes object. The axes object with title Numerical Solution contains an object of type patch.

Сравните этот результат с точным аналитическим решением и постройте ошибку.

p = model.Mesh.Nodes;
exact = (1 - p(1,:).^2 - p(2,:).^2)/4;
pdeplot(model,'XYData',u - exact')
title('Error');
xlabel('x')
ylabel('y')

Figure contains an axes object. The axes object with title Error contains an object of type patch.

Решения и ошибки с усовершенствованными сетками

Решите уравнение при совершенствовании mesh в каждой итерации и сравнении результата с точным решением. Каждое улучшение половины Hmax значение. Совершенствуйте mesh, пока норма по бесконечности вектора ошибок не будет меньше 5⋅10-7. 

hmax = 0.1;
error = [];
err = 1;
while err > 5e-7 % run until error <= 5e-7
    generateMesh(model,'Hmax',hmax); % refine mesh
    results = solvepde(model);
    u = results.NodalSolution;
    p = model.Mesh.Nodes;
    exact = (1 - p(1,:).^2 - p(2,:).^2)/4; 
    err = norm(u - exact',inf); % compare with exact solution
    error = [error err]; % keep history of err
    hmax = hmax/2;
end

Постройте норму по бесконечности вектора ошибок для каждой итерации. Значение ошибки уменьшается в каждой итерации.

plot(error,'rx','MarkerSize',12);
ax = gca;
ax.XTick = 1:numel(error);
title('Error History');
xlabel('Iteration');
ylabel('Norm of Error');

Figure contains an axes object. The axes object with title Error History contains an object of type line.

Постройте итоговую mesh и ее соответствующее решение.

figure
pdemesh(model); 
axis equal

Figure contains an axes object. The axes object contains 2 objects of type line.

figure
pdeplot(model,'XYData',u)
title('Numerical Solution');
xlabel('x')
ylabel('y')

Figure contains an axes object. The axes object with title Numerical Solution contains an object of type patch.

Сравните результат с точным аналитическим решением и постройте ошибку.

p = model.Mesh.Nodes;
exact = (1 - p(1,:).^2 - p(2,:).^2)/4;
pdeplot(model,'XYData',u - exact')
title('Error');
xlabel('x')
ylabel('y')

Figure contains an axes object. The axes object with title Error contains an object of type patch.

Документация Partial Differential Equation Toolbox

Поддержка