Note: This page has been translated by MathWorks. Click here to see

To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

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 example uses the `solvepde`

function. For the PDE Modeler app solution, see Poisson's Equation on Unit Disk: PDE Modeler App. Because the app and the programmatic workflow use different meshers, they yield slightly different results.

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