Exponenta Banner
exponenta event banner

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

Скрипт Open Live Script

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

Пример использует функцию solvepde. Для прикладного решения PDE Modeler смотрите уравнение Пуассона на Единичном диске: Приложение PDE Modeler. Поскольку приложение и программируемый рабочий процесс используют различный meshers, они приводят к немного отличающимся результатам.

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

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
Поддержка