Poisson's Equation on Unit Disk

This example shows how to numerically solve a Poisson's equation, compare the numerical solution with the exact solution, and refine the mesh until the solutions are close.

The Poisson equation on a unit disk with zero Dirichlet boundary condition can be written as $$-\Delta u=1$$ in $$\Omega$$, $$u=0$$ on $$\delta \Omega$$, where $$\Omega$$ is the unit disk. The exact solution is

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

For most PDEs, the exact solution is not known. However, the Poisson's equation on a unit disk has a known, exact solution that you can use to see how the error decreases as you refine the mesh.

Problem Definition

Create the PDE model and include the geometry.

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

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

figure 
pdegplot(model,"EdgeLabels","on"); 
axis equal

Specify zero Dirichlet boundary conditions on all edges.

applyBoundaryCondition(model,"dirichlet", ...
                             "Edge",1:model.Geometry.NumEdges, ...
                             "u",0);

Specify the coefficients.

specifyCoefficients(model,"m",0,"d",0,"c",1,"a",0,"f",1);

Solution and Error with a Coarse Mesh

Create a mesh with target maximum element size 0.1.

hmax = 0.1;
generateMesh(model,"Hmax",hmax);
figure
pdemesh(model); 
axis equal

Solve the PDE and plot the solution.

results = solvepde(model);
u = results.NodalSolution;
pdeplot(model,"XYData",u)
title("Numerical Solution");
xlabel("x")
ylabel("y")

Compare this result with the exact analytical solution and plot the error.

p = model.Mesh.Nodes;
exact = (1 - p(1,:).^2 - p(2,:).^2)/4;
pdeplot(model,"XYData",u - exact')
title("Error");
xlabel("x")
ylabel("y")

Solutions and Errors with Refined Meshes

Solve the equation while refining the mesh in each iteration and comparing the result with the exact solution. Each refinement halves the Hmax value. Refine the mesh until the infinity norm of the error vector is less than ${5\cdot 10}^{-7}$.

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 the infinity norm of the error vector for each iteration. The value of the error decreases in each iteration.

plot(error,"rx","MarkerSize",12);
ax = gca;
ax.XTick = 1:numel(error);
title("Error History");
xlabel("Iteration");
ylabel("Norm of Error");

Plot the final mesh and its corresponding solution.

figure
pdemesh(model); 
axis equal

figure
pdeplot(model,"XYData",u)
title("Numerical Solution");
xlabel("x")
ylabel("y")

Compare the result with the exact analytical solution and plot the error.

p = model.Mesh.Nodes;
exact = (1 - p(1,:).^2 - p(2,:).^2)/4;
pdeplot(model,"XYData",u - exact')
title("Error");
xlabel("x")
ylabel("y")