Home > Numerical Analysis, Partial Differential Equations > Finite Difference Method for 2D Laplace equation

Finite Difference Method for 2D Laplace equation

[Edit: This is, in fact Poisson’s equation.]

[For solving this equation on an arbitrary region using the finite difference method, take a look at this post.]

I will present here how to solve the Laplace equation using finite differences 2-dimensional case:

\displaystyle \begin{cases} - \Delta u=1 & \text{ in }\Omega=[0,1]^2 \\ u=0 &\text{ on }\partial \Omega. \end{cases}

Pick a step {h=1/N}, where {N} is a positive integer. The discretization of our function is a sequence of elements {u_{i,j}} with {0\leq i,j \leq N}. Then we use same the approximation for the derivative for both partial derivatives so that we have the approximation

\displaystyle -\Delta u(i/n,j/n)\simeq \frac{1}{h^2}(4u_{i,j}-u_{i,j-1}-u_{i,j+1}-u_{i+1,j}-u_{i-1,j}),

on the points situated in the interior of {\Omega} (namely for {1\leq i,j\leq N-1}).

The boundary conditions translate to {u_{0,j}=u_{N,j}=u_{i,0}=u_{i,N}=0} for {0\leq i,j \leq N}. Now on the points situated in the interior of {(0,1)} we have the relations:

\displaystyle \frac{1}{h^2}(4u_{i,j}-u_{i,j-1}-u_{i,j+1}-u_{i+1,j}-u_{i-1,j})=1,\ i=1..N-1,

and these translate to a system of equations

\displaystyle Au=f

where {A} is an {(N-1)^2 \times (N-1)^2} tridiagonal block matrix

\displaystyle A=\frac{1}{h^2}\text{tridiag}(-I,T,-I)


\displaystyle T=\begin{pmatrix} 4 & -1 & 0 &\cdots & 0 \\ -1 & 4 & -1 & \cdots & 0 \\ 0& -1 & 4 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 4 \end{pmatrix}=\frac{1}{h^2}\text{tridiag}(-1,4,-1)

and {f} is a column vector of length {(N-1)^2} of ones. Since {A} is invertible the presented system has a unique solution, which will be the desired approximation by finite differences of our function.

Below I present a simple Matlab code which solves the initial problem using the finite difference method. I also present a numerical example.

function laplacefd2(n);
U=reshape(u,[n-1 n-1])';
U=[zeros(n+1,1) [zeros(1,n-1); U; zeros(1,n-1)] zeros(n+1,1)];
[X,Y]= meshgrid(x,y); surf(X,Y,U);colorbar; colormap('cool');shading interp;

The function createmat assembles the matrix of the system which in the two dimensional case is a bit more complicated. Note that all the matrices are defined as sparse, i.e. with many zeros, and this helps a lot in the computation.

function A = createmat(n);
cell1 = sparse(diag(4*ones(n,1))+diag(-1*ones(n-1,1),-1)+diag(-1*ones(n-1,1),1));
cell2 =sparse(diag(-1*ones(n^2-n,1),-n)+diag(-1*ones(n^2-n,1),n));
A=sparse(tridiag + cell2);

  1. Dany
    December 26, 2012 at 12:48 pm

    Hi Beni,
    my question is how can i make constraints on some sells in my domain so that their value will allway stay “0”.
    lets say i have a 10×10 domain with 2 sources on both sides and in the middle i have a 2×2 block that the distribution should go around it.

    thank you


    • December 27, 2012 at 12:44 am

      I guess you could do that. To use the finite difference method the boundaries should be parallel to the coordonate axis. Putting a rectangular hole in the domain shouldn’t cause any problems.

      • Dany
        December 27, 2012 at 1:48 pm

        is there a different method where the boundries doesnt have to be parallel to the axis?
        like central difference method or maybe something else?

        thank you

      • Carla
        April 10, 2016 at 12:07 pm

        Hello 🙂
        I have a domain with a rectangular hole in the middle. However, I don’t know how to do the finite difference method with this. Could you give me some tips?
        Thank you 🙂

      • April 10, 2016 at 12:21 pm

        Hello. If you only need Dirichlet boundary conditions, you could use a penalization on the missing rectangle. The main ideas are presented here: https://mathproblems123.wordpress.com/2014/10/14/solving-poissons-equation-on-a-general-shape-using-finite-differences/

        If you do not have constraints on the software you want to use, you could try FreeFem, which can solve a variety of pdes on all kinds of domains.

  2. December 27, 2012 at 2:07 pm

    A method that works for domains of arbitrary shapes is the Finite Elements Method. The method is simple to describe, but a bit hard to implement. The Matlab PDE toolbox uses that method.

    The idea in 2D is: first discretize the boundary so that you have a polygonal shape. Choose a triangulation of this domain and consider only the values of your function in the triangle vertices. Then approximate the function by the piecewise affine function which is affine on each triangle and has the prescribed values in the vertices of the triangles. Then you write the equation in the weak formulation (using Sobolev spaces) and discretize your Sobolev-like space so that the differential equation turns to a linear system of equations.

    Maybe my presentation is a bit ambiguous, but it’s the best that I could do in a comment. Searching a bit you will find more details about the subject.

    The hardest part in the implementation is to use an algorithm to generate the triangulation and to turn your equation into a linear one.

    The finite difference method doesn’t work for domains with boundary parts non-parallel to the coordinate axis because it uses evenly spaced approximations for the partial derivatives in every direction. This will give distorted solutions near the parts of the boundary which are not parallel to the coordinate axis. It is also hard to write the linear system associated to a finite difference equation if your domain is not rectangular.

  3. mema
    September 1, 2013 at 9:43 pm

    i have un equally hy ,hx (grid mesh )how i program it with finit difference method in matlab

    ihave 110 by 32 matrixe

    thank you

    • February 23, 2014 at 9:12 pm

      You define an ordering on the matrix of points you have and this gives you a vector. Now you just need to find the matrix which multiplied with that vector produces the finite diferences. Maybe I will come to that in a following post.

  4. June 25, 2014 at 11:30 am

    hi …if i have to solve laplace’s equation for a square width and length = 2pi
    but the boundary conditions are (cos(x))^2 for the horizontal and (sin(x))^2 for vertical …how i can do it ?? and thank you 🙂

  5. October 14, 2014 at 2:40 pm

    How can I do this for an arbitary 2D domain ?

    • October 14, 2014 at 2:59 pm

      As you can read in come of my comments above, the finite difference method works well in square or rectangular domains. If your domain is arbitrary, the finite element method works. Matlab PDE tool uses that method.

      If you really want to use the finite difference method, you could penalize the points which are outside the domain by considering an operator of the form (-\Delta + \mu I) where \mu is zero inside your domain, and something very large outside. Then your matrix equation will transform to [(A+\text{diag}(1-\chi_\Omega)I]u=1.

    • October 14, 2014 at 4:31 pm
  6. October 29, 2014 at 2:30 pm

    Thanks man, I am looking for a way to solve helmholtz equation in 2D using this method.

  1. October 14, 2014 at 4:21 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 )

Twitter picture

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

Facebook photo

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

Google+ photo

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

Connecting to %s

%d bloggers like this: