Home > Numerical Analysis > Finite difference Method for 1D Laplace Equation

Finite difference Method for 1D Laplace Equation


I will present here how to solve the Laplace equation using finite differences. First let’s look at the one dimensional case:

\displaystyle \begin{cases} -u''=1 & \text{ in }[0,1] \\ u(0)=u(1)=0. \end{cases}

In this case it is easy to solve explicitly the differential equation and get that {\displaystyle u(x)=\frac{1-x^2}{2}}. However, sometimes it is not so easy or not possible at all to find an explicit solution, and there are some numerical methods which can be used to approximate our solution. In the one dimensional case, or in the case of a rectangular domain or a box in superior dimensions, the method of finite differences can be sucessfully used.

Pick a step {h=1/N}, where {N} is a positive integer. Then we use the approximations {u'(x)\simeq \displaystyle \frac{u(x+h)-u(x)}{h}} and {\displaystyle u''(x) \simeq \frac{u(x+h)+u(x-h)-2u(x)}{h^2}}.

We replace our function {u} by a discretization {u_k=u(k/N)} and the derivative {u''} becomes {-2u_i+u_{i+1}+u_{i-1}}. The boundary conditions translate to {u_0=u_N=0}. Now on the points situated in the interior of {(0,1)} we have the relations:

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

and these translate to a system of equations

\displaystyle Au=f

where {A} is an {N-1 \times N-1} matrix

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

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

We can change the boundary conditions to {u(0)=a} and {u(1)=b} by changing {f_1} and {f_{N-1}} using the formulas

\displaystyle f_1=1+N^2a,\ f_{N-1}=1+N^2 b.

In a similar way we can solve numerically the equation

\displaystyle \begin{cases} -u''=f & \text{ in }[0,1] \\ u(0)=u(1)=0. \end{cases}

where {f} is, for example, an arbitrary continuous function.

Below I present a simple Matlab code which solves the initial problem using the finite difference method and a few results obtained with the code.

function u = laplacefd1(n);
x=linspace(0,1,n+1);
A=sparse(diag(2*ones(n-1,1))+diag((-1)*ones(n-2,1),1)+diag((-1)*ones(n-2,1),-1));
left=0;
right=0;
f=ones(n-1,1);
f(1)=f(1)+(n)^2*left;
f(n-1)=f(n-1)+(n)^2*right;
u=(n)^2*A\f;

Different boundary conditions

Dirchlet boundary conditions

 

Advertisements
  1. Hannie Grace Jalon
    October 7, 2016 at 11:24 pm

    where do these formulas: \displaystyle f_1=1+N^2a,\ f_{N-1}=1+N^2 b come from??
    is this based from solving \displaystyle \frac{1}{h^2}(2u_i-u_{i+1}-u_{i-1})=1, for i=1 and i=N ??
    how did you get a “+N^2” in each formula?

    • October 8, 2016 at 1:07 am

      Yes, the formulas come from the discrete Laplacian formulation. Note that h=1/N so 1/h^2 = N^2. That’s why you get an N squared in the formulas.

      To see more precisely why it is enough to change f_1,f_N just write the discrete Laplacian for u_0,u_1,u_2. You’ll get N^2(2u_1-u_0-u_2) = f_1. Now if you pass the term with u_0 to the other side you see that you have the same matrix equation like in the zero boundary condition, but with different f.

  2. David
    February 8, 2017 at 2:28 pm

    Hi, in your first equation you have the problem -u”=1 with boundary conditions u(0) = u(1) = 0. You state that the analytical solution to this problem is 0.5*(1-x^2).

    Now when I insert the boundary condition into your analytical solution, I get u(0) = 0.5*(1-0^2) = 0.5.

    Are you sure the boundary conditions in your problem are correct?

    • February 22, 2017 at 4:51 pm

      You are right. The solution in my post is, in fact, on the interval [-1,1]. In order for it to work on [0,1] you need to change it a bit. I guess 0.5*x*(x-1) should work instead.

  1. No trackbacks yet.

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: