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