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