В этом примере показано, как численно решить уравнение Пуассона, сравните числовое решение с точным решением и совершенствуйте mesh, пока решения не близки.
Уравнение Poisson на единичном диске с нулем граничное условие Дирихле может быть записано как \in , on , где единичный диск. Точное решение
Для большинства УЧП не известно точное решение. Однако уравнение Пуассона на единичном диске имеет известное, точное решение, которое можно использовать, чтобы видеть, как ошибка уменьшается, когда вы совершенствовали 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, пока норма по бесконечности вектора ошибок не будет меньше .
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')
