This is machine translation

Translated by Microsoft
Mouseover text to see original. Click the button below to return to the English version of the page.

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.

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 Solve Poisson's Equation on a Unit Disk: PDE Modeler App. The problem formulation is in , on , where is the unit disk. The exact solution is

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

Specify zero Dirichlet boundary conditions on all edges.


Specify the coefficients.


Create a mesh with target maximum element size 0.1.


Solve the PDE and plot the solution:

results = solvepde(model);
u = results.NodalSolution;
title('Numerical Solution');

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

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;

The value of the error decreases in each iteration.

ax = gca;
ax.XTick = 1:numel(error);
title('Error History');
ylabel('Norm of Error');

Was this topic helpful?