Documentation Center

  • Trial Software
  • Product Updates

Poisson's Equation with Point Source and Adaptive Mesh Refinement

This example solves a Poisson's equation with a delta-function point source on the unit disk using the adaptmesh function in the Partial Differential Equation Toolbox™.

Specifically, we solve the Poisson's equation

$$-\Delta u=\delta(x,y)$$

on the unit disk with zero Dirichlet boundary conditions. The exact solution expressed in polar coordinates is

$$u(r,\theta)=\frac{\log(r)}{2\pi},$$

which is singular at the origin.

By using adaptive mesh refinement, the Partial Equation Toolbox™ can accurately find the solution everywhere away from the origin.

Problem Definition

The following variables will define our problem:

  • g: A specification function that is used by initmesh. For more information, please see the documentation page for circleg and pdegeom.

  • b: A boundary file used by assempde. For more information on this file, please see the documentation pages for circleb1 and pdebound.

  • c, a: The coefficients of the PDE.

  • f: A function that captures a point source at the origin. It returns 1/area for the triangle containing the origin and 0 for other triangles.

g='circleg';
b='circleb1';
c=1;
a=0;
f='circlef';

Generate Mesh

adaptmesh solves elliptic PDEs using adaptive mesh generation. The 'tripick' parameter lets the user specify a function that returns which triangles will be refined in the next iteration. circlepick returns triangles with computed error estimates greater a given tolerance. The tolerance is provided to circlepick using the 'par' parameter.

[u,p,e,t]=adaptmesh(g,b,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.

Plot Finest Mesh

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

Plot Error

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');

Plot FEM Solution on Finest Mesh

figure;
pdeplot(p,e,t,'xydata',u,'zdata',u,'mesh','off');

Was this topic helpful?