MATLAB Examples

Solve Poisson's Equation on a Unit Disk

This example shows how to solve Poisson's equation using the programmatic workflow. For the PDE Modeler app solution, see docid:pde_ug.bvhf75n. The problem formulation is $- \Delta u = 1$ in $\Omega$, $u = 0$ on $\delta \Omega$, where $\Omega$ is the unit disk. The exact solution is

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

The code compares the solution obtained with the PDE Toolbox with the exact solution, and refines the mesh until the solutions are close.

First, create a function that parameterizes the 2-D geometry — in this case, a unit circle. The circleg.m file contains a function that returns the coordinates of points on the unit circle's boundary. You can display the file by using the command type circleg.

Create the PDE model and include the geometry.

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

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

Create a mesh with target maximum element size 0.1.

generateMesh(model,'Hmax',0.1);

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

Now use the coarser mesh, refining it in each iteration. Solve the equation using the original coarse and refined meshes, each time comparing the result with the exact solution.

hmax = 1;
error = [];
err = 1;
while err > 0.001 % run until error <= 0.001
    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

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