Partial Differential Equation Toolbox

Poisson's Equation on a Unit Disk

This example shows how to numerically solve a Poisson's equation using the assempde function in Partial Differential Equation Toolbox™.

The particular PDE is

$$ -\Delta u = 1$$

on the unit disk with zero-Dirichlet boundary conditions. The exact solution is

$$ u(x,y) = \frac{1-x^2-y^2}{4}. $$

For most partial differential equations, the exact solution is not known. In this example, however, we can use the known, exact solution to show how the error decreases as the mesh is refined.

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 pages for circleg and pdegeom.

  • c, a, f: The coefficients and inhomogeneous term.

g = @circleg;
c = 1;
a = 0;
f = 1;

Boundary Conditions

Plot the geometry and display the edge labels for use in the boundary condition definition.

figure
pdegplot(g, 'edgeLabels', 'on');
axis equal
% Create a pde entity for a PDE with a single dependent variable
numberOfPDE = 1;
pb = pde(numberOfPDE);
% Create a geometry entity
pg = pdeGeometryFromEdges(g);
% Solution is zero at all four outer edges of the circle
pb.BoundaryConditions = pdeBoundaryConditions(pg.Edges(1:4), 'u', 0);

Generate Initial Mesh

The function initmesh takes a geometry specification function and returns a discretization of that domain. The 'hmax' option lets the user specify the maximum edge length. In this case, because the domain is a unit disk, a maximum edge length of one creates a coarse discretization.

[p,e,t] = initmesh(g,'hmax',1);
figure
pdemesh(p,e,t);
axis equal

Refinement

We repeatedly refine the mesh until the infinity-norm of the error vector is less than a $10^{-3}$.

For this domain, each refinement halves the lengths of the edges of the triangles that compose the mesh. Note also that the error decreases by a factor of about one-third.

er = Inf;
while er > 0.001
    [p,e,t] = refinemesh(g,p,e,t);
    u = assempde(pb,p,e,t,c,a,f);
    exact = (1-p(1,:).^2-p(2,:).^2)'/4;
    er = norm(u-exact,'inf');
    fprintf('Error: %e. Number of nodes: %d\n',er,size(p,2));
end
Error: 1.292265e-02. Number of nodes: 25
Error: 4.079923e-03. Number of nodes: 81
Error: 1.221020e-03. Number of nodes: 289
Error: 3.547924e-04. Number of nodes: 1089

Plot Final Mesh

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

Plot Error

The scale of the vertical axis shows that the error is small and of the order $10^{-4}$.

figure
pdesurf(p,t,u-exact);

Plot FEM Solution

figure
pdesurf(p,t,u);