Уравнение Пуассона с точечным источником и адаптивное улучшение Mesh

В этом примере показано, как решить уравнение Пуассона с источником единицы функциональности дельты на единичном диске с помощью adaptmesh функция.

А именно, решите уравнение Пуассона

на единичном диске с нулем граничные условия Дирихле. Точное решение, описанное в полярных координатах,

который сингулярен в начале координат.

При помощи адаптивного улучшения mesh Уравнение в частных производных Toolbox™ может точно найти решение везде далеко от источника.

Следующие переменные описывают задачу:

c = 1;
a = 0;
f = @circlef;

Создайте Модель УЧП с одной зависимой переменной.

numberOfPDE = 1;
model = createpde(numberOfPDE);

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

g = @circleg;
geometryFromEdges(model,g);

Постройте геометрию и отобразите метки ребра.

figure; 
pdegplot(model,'EdgeLabels','on'); 
axis equal
title 'Geometry With Edge Labels Displayed';

Задайте нулевое решение во всех четырех внешних краях круга.

applyBoundaryCondition(model,'dirichlet','Edge',(1:4),'u',0);

adaptmesh решает эллиптические УЧП с помощью адаптивной генерации mesh. tripick параметр позволяет вам задать функцию, которая возвращается, какие треугольники будут усовершенствованы в следующей итерации. circlepick возвращается треугольники с вычисленной ошибкой оценивает больше данный допуск. Допуск предоставляется circlepick использование параметра 'паритета'.

[u,p,e,t] = adaptmesh(g,model,c,a,f,'tripick', ...
                                    'circlepick', ...
                                    'maxt',2000, ...
                                    'par',1e-3);

Number of triangles: 258
Number of triangles: 515
Number of triangles: 747
Number of triangles: 1003
Number of triangles: 1243
Number of triangles: 1481
Number of triangles: 1705
Number of triangles: 1943
Number of triangles: 2155

Maximum number of triangles obtained.

Постройте самую прекрасную mesh.

figure; 
pdemesh(p,e,t); 
axis equal

Постройте ошибочные значения.

x = p(1,:)';
y = p(2,:)';
r = sqrt(x.^2+y.^2);
uu = -log(r)/2/pi;
figure;
pdeplot(p,e,t,'XYData',u-uu,'ZData',u-uu,'Mesh','off');

Постройте решение FEM на самой прекрасной mesh.

figure;
pdeplot(p,e,t,'XYData',u,'ZData',u,'Mesh','off');