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.
Figure 1. Problem definition (left) and domain discretization (right)
As an example, let’s consider the Poisson equation, . 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 . From Taylor’s Series, we know that a value of a function some distance from a known point x can be estimated from derivatives as
In other words, the second derivative at x given by
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
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:
The first and third derivative conveniently cancel out, and the resulting expression is second-order accurate.
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.
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 or . From the Taylor expansion, the first derivative is given by
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
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.