## 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 which contains the shape and add a penalization on the exterior of your domain . The problem to solve becomes something like:

Note that doing this we do not need to impose the boundary condition on . This is already imposed by , and the fact that is forced to be zero outside .

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

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

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

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

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 will be given by its radial function written in a finite Fourier series.

The program takes as parameters.

function penalized_laplace(phi,n,M); form = phi(:); col = zeros(n^2,1); col(1)=4; col(2)=-1; col(n+1)=-1; lin = col'; A = n^2/9*sptoeplitz(col,lin); phi = phi_fourier(form,[n n]); phi = 1-phi(:); G=spdiags(phi,0,n^2,n^2); 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 = []; end x=linspace(-1.5,1.5,nx(1)); y=linspace(-1.5,1.5,nx(2)); 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)'; temp=zeros(nx(1),nx(2)); 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)) temp(i,j)=1; end; end; end; res=temp;

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.

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

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