Exponenta Banner
exponenta event banner

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

Открыть сценарий в реальном времени

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

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

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

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

Определение проблемы

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

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);

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

Создайте сетку с целевым максимальным размером элемента 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.

Решите PDE и постройте график решения.

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.

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

Решите уравнение, уточняя сетку в каждой итерации и сравнивая результат с точным решением. Каждое уточнение наполовину Hmax значение. Уточняйте сетку до тех пор, пока бесконечная норма вектора ошибок не станет меньше 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.

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

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.

Документация по набору инструментов для дифференциальных уравнений в частных производных

Поддержка