Home > Numerical Analysis, Partial Differential Equations, Sobolev Spaces > Solving Poisson’s equation on a general shape using finite differences

Solving Poisson’s equation on a general shape using finite differences

One of the questions I received in the comments to my old post on solving Laplace equation (in fact this is Poisson’s equation) using finite differences was how to apply this procedure on arbitrary domains. It is not possible to use this method directly on general domains. The main problem is the fact that, unlike in the case of a square of rectangular domain, when we have a general shape, the boudary can have any orientation, not only the orientation of the coordinate axes. One way to avoid approach this problem would be using the Finite Element Method. Briefly, you discretize the boundary, you consider a triangulation of the domain with respect to this discretization, then you consider functions which are polynomial and have support in a few number of triangles. Thus the problem is reduced to a finite dimensional one, which can be written as a matrix problem. The implementation is not straightforward, since you need to conceive algorithms for doing the discretization and triangulation of your domain.

One other approach is to consider a rectangle {D} which contains the shape {\Omega} and add a penalization on the exterior of your domain {\Omega}. The problem to solve becomes something like:

\displaystyle (-\Delta +\mu I) u = 1 \text{ on }D

where {\mu} is defined by

\displaystyle \mu(x) = \begin{cases} 0 & x \in \Omega \\ + \infty & x \notin \Omega\end{cases}. \ \ \ \ \ (1)

Note that doing this we do not need to impose the boundary condition {u=0} on {\partial \Omega}. This is already imposed by {\mu}, and the fact that {u} is forced to be zero outside {\Omega}.

The above consideration can be made rigorous, if we define the relaxation of a Dirichlet problem. For {\mu} a capacitary measure, we can consider the problem

\displaystyle \begin{cases} -\Delta u + u \mu = f,\\ u \in H_0^1(D) \cap L^2(\mu).\end{cases}

Using classical arguments like the Lax-Milgram theorem, it can be proved that the above problem has a solution. Moreover, in the case where {\mu} has the form (1), it is possible to prove that the solution is in {H_0^1(\Omega)}.

This relaxation is useful if we want to compute numerically the solution of the Poisson equation

\displaystyle -\Delta u = 1,\ u \in H_0^1(\Omega)

using finite differences. As stated, we consider the shape included in a square and define the measure {\mu} to be zero inside {\Omega} and a very big number outside. The discretized problem becomes

\displaystyle (A+\text{diag}(1-\chi_\Omega)M \cdot I)u = 1.

This is just a system of equations, which can be solved right away.

I write below the Matlab code to solve such a problem in two dimensions. The shape {\Omega} will be given by its radial function written in a finite Fourier series.

\displaystyle \rho(\theta) = a_0 +\sum_{k=1}^N a_k \cos(kx)+b_k\sin(kx).

The program takes {(a_0,a_1,...,a_N,b_1,...,b_N)} as parameters.

 function penalized_laplace(phi,n,M);
form = phi(:);
col = zeros(n^2,1);
lin = col';
A = n^2/9*sptoeplitz(col,lin);
phi = phi_fourier(form,[n n]);
phi = 1-phi(:);
u = (A+M*G)\ones(n^2,1);
[X,Y] = meshgrid(linspace(-1.5,1.5,n),linspace(-1.5,1.5,n));
surf(X,Y,reshape(u,[n n]));%view(2);
shading interp;
axis([-1.5 1.5 -1.5 1.5 0 max(u)]) 

The function phi_fourier.m computes the characteristic function on the grid starting from the Fourier series.

function [res]=phi_fourier(vec,nx,opts);

if nargin < 3
opts = [];

vec = vec(:);

dim = length(vec);
m = (dim-1)/2;

a0 = vec(1);
a = vec(2:m+1);
b = vec(m+2:2*m+1);

ind = (1:m)';
for i= 1:nx(1)
for j= 1:nx(2)
theta = cart2pol(x(i),y(j));
if norm([x(i); y(j)])<=a0 + dot(a,cos(theta*ind))+dot(b,sin(theta*ind))

Here are some of the results you can obtain.

penalized_laplace([1 0 0 0 0],100,1e6) solves the Poisson equation on the unit disk.


penalized_laplace([1 0 0 0 0.3],100,1e6); view(2) produces the result presented in the second figure.


  1. October 29, 2014 at 2:25 pm

    Hi man, Thanks for the code. After downloading sptoeplitz function from the file exchange the code worked. I assume that this method is applicable to laplace eigenvalues also (2D Helmholtz equation).

    • October 29, 2014 at 6:53 pm

      I’m glad you managed to make the code work. Yes, the method also works for the eigenvalues of the Laplace operator.

  1. October 27, 2014 at 7:04 pm

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s

%d bloggers like this: