## 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:

Pick a step , where is a positive integer. The discretization of our function is a sequence of elements with . Then we use same the approximation for the derivative for both partial derivatives so that we have the approximation

on the points situated in the interior of (namely for ).

The boundary conditions translate to for . Now on the points situated in the interior of we have the relations:

and these translate to a system of equations

where is an tridiagonal block matrix

where

and is a column vector of length of ones. Since 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);`

x=linspace(0,1,n+1);

y=linspace(0,1,n+1);

A=createmat(n-1);

f=ones((n-1)^2,1);

u=n^2*A\f;

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

tridiag=kron(eye(n),cell1);

cell2 =sparse(diag(-1*ones(n^2-n,1),-n)+diag(-1*ones(n^2-n,1),n));

A=sparse(tridiag + cell2);

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

Dany

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.

ok.

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

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 🙂

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.

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.

hi

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

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.

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 🙂

You need to just change the boundary values. I’ll try and do it when I have some spare time.

Hi …ok i appreciate it ..but may i ask if you could do it as soon as possible ? Because i have a university project to present before the end of the month ..and thank you …again …

It is not my style to do someone else’s work, plus, I do not have the time now… If it is your project, it should be your work. Take a look at the document in the next link: http://www.pa.uky.edu/~plaster/phy611/Laplace_Finite_Difference.pdf or at this program: http://www.mathworks.ch/matlabcentral/fileexchange/38091-2d-laplace-equation/content/Laplace_equation_2D.m

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

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 where is zero inside your domain, and something very large outside. Then your matrix equation will transform to .

Meanwhile I’ve written a post where I talk about this. Here it is: https://mathproblems123.wordpress.com/2014/10/14/solving-poissons-equation-on-a-general-shape-using-finite-differences/

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