Main Content

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.

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

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

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