Exponenta Banner

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

Открыть Live Script

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

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

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

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

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

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

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

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

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

Figure contains an axes. The axes 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

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

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

Figure contains an axes. The axes 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. The axes 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. The axes 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. The axes with title Error History contains an object of type line.

Постройте график конечного mesh и ее соответствующего решения.

figure
pdemesh(model); 
axis equal

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

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

Figure contains an axes. The axes 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. The axes with title Error contains an object of type patch.

Partial Differential Equation Toolbox производными

Поддержка