PIC-C, Particle In Cell Consulting, LLC. Plasma physics and rarefied gas simulation codes.

The Finite Difference Method

Numerical simulations of physical processes generally involve solving some differential equation on a computational domain too complicated to solve analytically. Solving simple systems by “hand” is quite possible in one-dimension. But things get more complicated as you go to higher dimensions. If the domain has a nice shape – that is, if it’s rectangular or cylindrical in nature – you may be able to solve it analytically using techniques such as separation of variables. But as you start introducing irregularities in the boundary or in the forcing function, things start getting hairy really soon. In that case, going to a numerical solution is the only viable option.

The Finite Difference Method (FDM) is a way to solve differential equations numerically. It is not the only option, alternatives include the finite volume and finite element methods, and also various mesh-free approaches. However, FDM is very popular. The popularity of FDM stems from the fact it is very simple to both derive and implement numerically.

finite difference simulation domain finite difference discretization

Figure 1. Problem definition (left) and domain discretization (right)

As an example, let’s consider the Poisson equation,  \nabla^2\phi=-\dfrac{\rho}{\epsilon_0}. This equations governs the variation of electric potential given some charge density distribution. It is one of the most fundamental equations in the field of electro-static plasma simulations. We want to solve this equation numerically on a rectangular domain shown in Figure 1 subject to the boundaries listed in the figure. The domain contains of two regions of fixed potential along upper and bottom edge – these could represent charged electrodes. Remaining edges have zero electric field, except for the left edge, on which electric field is specified.

We start by discretizing the domain – in other words, overlaying a computational mesh over the domain. In the Finite Difference method, solution to the system is known only on on the nodes of the computational mesh. As such, it is important to chose mesh spacing fine enough to resolve the details of interest. In addition, cell edges must coincide with the axis of the coordinate system being used. This is one of the main disadvantages of FDM: complex geometries cannot be directly resolved by fitting the mesh to the object boundary.

Finite Difference representation of derivatives

We are looking for the solution to  \nabla^2\phi=\dfrac{\partial^2\phi}{\partial x^2}+\dfrac{\partial^2\phi}{\partial y^2}=-\dfrac{\rho}{\epsilon_0}. From Taylor’s Series, we know that a value of a function some distance \Delta x from a known point x can be estimated from derivatives as
f(x+\Delta x)=f(x)+\dfrac{f'(x)}{1!}(\Delta x) + \dfrac{f''(x)}{2!}(\Delta x)^2 + \dfrac{f'''(x)}{3!}(\Delta x)^3+ O(4)+\ldots
In other words, the second derivative at x given by
f''(x) = \dfrac{2}{\Delta^2x}\left[ f(x+\Delta x) - f(x) - f'(x)\Delta x -  \dfrac{f'''(x)}{6}\Delta^3 x+ O(4)+\ldots \right]
The expression above is known as the forward difference. We are estimating the derivative at a point using data to the front (in positive direction) of that point. We can obtain a similar expression by going backward
f(x-\Delta x)=f(x)+\dfrac{f'(x)}{1!}(-\Delta x) + \dfrac{f''(x)}{2!}(-\Delta x)^2 + \dfrac{f'''(x)}{3!}(-\Delta x)^3+ O(4)+\ldots
f''(x) = \dfrac{2}{\Delta^2x}\left[ f(x-\Delta x) - f(x) + f'(x)\Delta x +  \dfrac{f'''(x)}{6}\Delta^3 x+ O(4)+\ldots \right]

We now have two expressions for the second derivative at point x. Instead of using one of the other, it’s best practice to use their average. By adding the two expression and dividing by two, we obtain the central difference representation for the second derivative:
f''(x) = \dfrac{f(x-\Delta x) - 2f(x) + f(x+\Delta x)}{\Delta^2x} + O(4) +\ldots
The first and third derivative conveniently cancel out, and the resulting expression is second-order accurate.

Numerical Implementation

On the computational domain, potential is no longer a continuous function, instead, it is given by a collection of values at node indices as phi[i][j]. Assuming uniform spacing between the grid nodes in each directions, the second derivative in the x direction is written as d2f_dx2 = (phi[i-1][j]-2*phi[i][j]+phi[i+1][j])/(dx*dx). Adding the expression for the y direction, we obtain the discretized version of the Poisson equation:
(phi[i-1][j]-2*phi[i][j]+phi[i+1][j])/(dx*dx) + (phi[i][j-1]-2*phi[i][j]+phi[i][j+1])/(dy*dy) = -rho[i][j]/eps0

That’s it, quite simple. The expression above gives us an expression that can be used to solve for potential everywhere inside the domain. But to complete the problem, we need to include the boundaries.

Boundary Conditions

Our problem has two types of boundary conditions: fixed potential along portions of top and bottom boundary, and fixed derivative (electric field) on the remaining nodes.

The first case is known as Dirichlet boundary condition. It is simple to implement. On each node along the boundary we have phi[i][j]=g[i][j]

The second condition is known as Neumann. On the left face, we have \nabla \phi \cdot (-i)=-E_0 or \dfrac{\partial \phi}{\partial y}=E_0. From the Taylor expansion, the first derivative is given by
f'(x)=\dfrac{f(x+\Delta x)-f(x)}{\Delta x} + O(2) + \ldots
This is the forward difference for the first derivative. Central difference can be derived using process analogous to above. However, since we do not have any data in the backward direction along the y=0 face, we are forced to use this less-accurate representation. This equation is implemented numerically as (phi[i+1][j]-phi[i][j])/dx = E0

Numerical implementation

This simple example implemented in a simple Java Finite Difference Solver. Give the program a try (I used Eclipse to write it) and let me know if you have any questions.

Figure 2. Solution to the sample problem using the Finite Difference Method.

For an alternative, see The Finite Volume Method.

Don’t forget to link to this article if you find it useful.

Subscribe to the newsletter and follow us on Twitter. Send us an email if you have any questions.

23 comments to “The Finite Difference Method”

  1. lubos
    March 1, 2011 at 6:59 pm

    Please leave a comment if you have any questions. Thanks!

    • Dimitris Leivaditis
      September 14, 2014 at 2:19 am


      I wonder how B[u]=1/(dx)^2, T[u]=1/(dy)^2 etc. were derived. I know it has something to do with the coefficient matrix but I cannot figure this out. For example, if the equation was again Poisson equation but the second derivative was multiplied by x^2, how do you compute these quantities?

      Thanks in advance,


      • Dimitris Leivaditis
        September 16, 2014 at 12:17 am

        Ok, I figured it out. They are just the coefficients of the top, bottom, left, right and central node after you discretize the equation.

        • September 18, 2014 at 9:24 pm

          Great, and sorry for not getting back to you in time.

  2. alex
    April 13, 2011 at 8:07 am

    Have you considered making your java code an applet so that it can be embedded into the website?

  3. John
    August 5, 2011 at 10:08 am

    Here is my web based javascript version of this solver.


    • August 5, 2011 at 10:41 am

      Wow, you rock! This is totally awesome. Here is a next task for you :) Use VTK to extract isosurface from some 3D field, save as VRML, and show as 3D object.

  4. August 6, 2011 at 6:56 am

    Here is my 3D version of solver.


    I haven’t had the chance to look at VTK to find out what they do with isosurface. All I did was plot the computational mesh points calculated by solver in 3D.

  5. john
    August 18, 2011 at 1:28 pm

    Can you verify the neumann boundary conditions in your solver.java .

    /* set neumann boundaries along the edges */
    for (i = 0; i < nx; i++) {
    /* bottom */
    T[i] = 1 / dx;
    C[i] = -1 / dx;

    /* top */
    B[(ny – 1) * nx + i] = 1 / dx;
    C[(ny – 1) * nx + i] = -1 / dx;

    Should the neumann boundary along the top and bottom be normal to the boundary. Does it make more sense to use 1 / dy instead of 1 / dx for the top and bottom boundary in the above equations. Also, are you supposed to perform a dot product with the (external) normal to the boundary, so does the top get multiplied by +1 (+y direction) and the bottom get multiplied -1 (-y direction) according to this wiki.


  6. john
    August 20, 2011 at 4:42 pm

    I think the equations in section 1.2 of the article dealing with insulated neumann boundaries du / dn = 0 is not quite right. I think the term in 1.2.2 for instance (hy / hx) * Uend-1,j should be added instead of subtracted. The same applies to the other insulated boundary equations in section 1.2 . Doing so gives a result that is similar to the one in this article. I also used the equations for the corner points from the boundary article for this latest result.


  7. john
    August 24, 2011 at 12:21 pm

    I modified the solver routine to use a Full Multigrid Cycle to solve the problem. Here is the link to the solution using multigrid.


  8. October 24, 2011 at 11:22 am

    Here is an example of the Finite Difference Time Domain method in 1D which makes use of the leapfrog staggered grid.


    You can learn more about the fdtd method here.


    and here.


  9. October 27, 2011 at 10:25 am

    Here is an example of the Finite Difference Time Domain method in 2D.


  10. derek
    November 16, 2011 at 2:46 am

    Hello again! You have many interesting things on this site. I have two questions please: 1) Is the value -rho[i][j]/eps0 equal to zero for every cell? 2) And what value of dx did you use to get your nice colorful graph? I used dx = 0.01 because anything bigger Excel doesn’t like it (or it seemed so…). Thanks a lot for your generous time!!!

    • November 16, 2011 at 9:13 am

      This depends on the local charge density. I believe it is zero in this particular example but in general it won’t be. Plasma is generally quasineutral, which means that rho is approximately zero at spatial scales greater than the Debye length. However, for a Poisson solver to converge, the cell size must be smaller or at least the order of the local Debye length, and at that spatial scale you will get some small deviations from quasineutrality. These deviations are the source of the restoring force that keeps ions and electrons close to each other.

      About your second question, these plots were actually generated by Tecplot. Tecplot is a commercial program for data visualization. It is a very nice package, but it seems that every year it picks up more features you don’t need and the price goes correspondingly up. There are some free data visualization alternatives available, such as Paraview and VisIt. I have also developed my own visualization tool many years back called capVTE. Unfortunately, it does not currently work anymore, but I am hoping to revive it shortly.

  11. John
    January 12, 2012 at 11:25 am
  12. John
    January 15, 2012 at 11:42 am

    I tweaked my matlab solution a little playing around with the boundary conditions. This is what matlab gives me when the left neumann boundary condition is 10. (First derivative in x direction set to 10)


    And this is what I get when I set the left neumann boundary condition to -2 .


  13. John
    May 8, 2012 at 9:34 am

    I tried implementing this example using the Finite Element Method in 2D on MATLAB. This is the result that was produced.


    I based it on the FEM2D code described is this paper.


    Finite Element MATLAB code can be found here.


  14. John
    May 11, 2012 at 10:38 am

    Here is another implementation using Finite Elements to solve in MATLAB and distmesh to generate the mesh.


    I used distmesh to automatically generate a mesh of triangles.


    I also modified this MIT MATLAB code for a 2D Finite Element solver to read in the mesh data, add Neumann boundary conditions and solve the problem.


  15. John
    May 11, 2012 at 11:05 am

    Here is one more example using 2D Finite Elements with distmesh on a disk shape object with boundary conditions of 0.5 and 7.0 on the inner and outer boundaries of the disk.


Leave a Reply

You can use these tags: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong> In addition, you can use$latex ...$ to include equations.