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

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

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

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

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

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

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

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

Решите уравнение при совершенствовании 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');

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

figure
pdemesh(model); 
axis equal

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

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

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

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

Поддержка