Finite differencing: Introduction

From Interactive System for Ice sheet Simulation
Jump to: navigation, search


(There is a PDF of the slides relating to this section available here, but this page contains all the material and more...)

Some notation

Finite differences work by dividing a continuous variable up into a succession of points, or nodes. For simplicity, we'll think about equally-spaced nodes in two spatial dimensions, and evolving in time. We can define a set of nodes in time and space:

\left\{(i\Delta x, j\Delta y, n\Delta t)\right\},

where i, j and n are the indices of the nodes, \Delta x and \Delta y are the distances between the nodes (the grid-spacing), and \Delta t is the timestep. A continuous function can be mapped onto the grid as follows:

\phi(x,y,t)\rightarrow\phi(i\Delta x, j\Delta y, n\Delta t).

For convenience, the discrete function can be noted using subscripts for the spatial indices, and a superscript for time[1]:

\phi(i\Delta x, j\Delta y, n\Delta t)\equiv\phi_{i,j}^n.

This may seem like a potential source of confusion, but the equations we're going to be solving don't generally involve higher powers of the quantities being differentiated. In practice, this compact notation is very handy!

My first finite difference expression

The nice thing about finite differences is that they're very intuitive. Consider a first-order spatial derivative:

\frac{\partial \phi}{\partial x}.

If we had a plot of \displaystyle\phi in front of us, and we wanted to calculate a rough-and-ready value of this derivative, we might read off the values of \displaystyle\phi at two nearby points (let's call them x_1 and x_2), and calculate the gradient from that:

\frac{\partial \phi}{\partial x}\approx\frac{\phi(x_2)-\phi(x_1)}{x_2-x_1}

Of course, these two points could be on the grid we defined earlier, and written using our compact notation:

\frac{\partial \phi}{\partial x}\approx\frac{\phi_{i+1}-\phi_i}{\Delta x}

So, there we have it: our first finite difference expression. But wait — there are some unanswered questions:

  • We've calculated a derivative, but where is it valid? Is this the slope at point x_i or point x_{i+1}? Or somewhere in between?
  • We know this is an approximation to the real value of the derivative, but how much of one? Can we quantify its accuracy?

To answer these questions, we need to look at finite differences in a bit more depth, using Taylor's Theorem

Deriving Finite Differences

Taylor's Theorem allows the approximation of a differentiable function at a given point by a polynomial whose coefficients depend on the function's derivatives at that point. For a differentiable function f, Taylor's Theorem is:

 f(x) = f(a) + \frac{(x - a)}{1!}\left.\frac{df}{dx}\right|_a + \frac{(x - a)^2}{2!}\left.\frac{d^2f}{dx^2}\right|_a + \cdots + \frac{(x - a)^n}{n!}\left.\frac{d^nf}{dx^n}\right|_a + R_n(x).

Here, a is the point where the derivatives are evaluated, and x is the point of interest. R_n(x) is the truncation error, since the series given here is of finite length. For most finite difference applications, second-order accuracy is sufficient. We can rewrite the series in terms of our grid, truncating the expression at the second-order term:

 \phi_{i+1} = \phi_i + \frac{\Delta x}{1!}\left.\frac{d\phi}{dx}\right|_i + \frac{\Delta x^2}{2!}\left.\frac{d^2\phi}{dx^2}\right|_i + \mathcal{O}(\Delta x^3) .

The truncation error has been rewritten to show that it is a third-order term. To get a second-order discretization of our spatial derivative, we also need the Taylor Series for \displaystyle\phi_{i-1}:

 \phi_{i-1} = \phi_i - \frac{\Delta x}{1!}\left.\frac{d\phi}{dx}\right|_i + \frac{\Delta x^2}{2!}\left.\frac{d^2\phi}{dx^2}\right|_i + \mathcal{O}(\Delta x^3) .

By subtracting the second equation from the first, we arrive at

 \phi_{i+1} - \phi_{i-1} = 2\Delta x \left.\frac{d\phi}{dx}\right|_i + \mathcal{O}(\Delta x^3),

which can be rearranged finally to give

 \left.\frac{d\phi}{dx}\right|_i = \frac{\phi_{i+1} - \phi_{i-1}}{2\Delta x} + \mathcal{O}(\Delta x^3).

So, for this finite difference expressions, we have answered both questions posed above using Taylor's Theorem: we know where this approximation is valid (at x_i) and how accurate it is (second-order accurate, with 'error' terms of order \Delta x ^3).

Stencils, and more derivatives

The expression we've just derived is a centered-difference expression, as the points it relies upon (x_{i+1} and x_{i-1}) lie either side of the point where we're evaluating the derivative (x_i). This set of points is known as the expression's stencil. It's important to realise that many different stencils can be used to calculate a given derivative, and each entails manipulating Taylor Series expressions in a different way. With this in mind, we return to the original finite difference expression we devised:

\frac{\partial \phi}{\partial x}\approx\frac{\phi_{i+1}-\phi_i}{\Delta x}.

As we noted above, there's no indication in the expression as to where the derivative is valid. It turns out that this is important with regard to its accuracy. First, let's consider this as an off-centred expression, evaluated at x_i:

\left.\frac{\partial \phi}{\partial x}\right|_i\approx\frac{\phi_{i+1}-\phi_i}{\Delta x}.

We can easily rearrange this into the form of a Taylor Series:

\phi_{i+1} \approx \phi_i + \frac{\Delta x}{1!}\left.\frac{\partial \phi}{\partial x}\right|_i.

By comparing with the Taylor Series expressions given earlier, it becomes clear that this has one fewer terms: the truncation error is second-order, and so the expression is of first-order accuracy:

\phi_{i+1} = \phi_i + \frac{\Delta x}{1!}\left.\frac{\partial \phi}{\partial x}\right|_i + \mathcal{O}(\Delta x^2).

However, we can also think of the derivative being evaluated halfway between x_i and x_{i+1}:

\left.\frac{\partial \phi}{\partial x}\right|_{i+\frac12}\approx\frac{\phi_{i+1}-\phi_i}{\Delta x}.

Of course, this is now a centered difference, and can be constructed in a similar manner to the centred difference expression given previously. For similar reasons, it is also second-order accurate.

Second derivatives and associated error terms are also computed from Taylor's series, this time considering sums of series

 \phi_{i+1} = \phi_i + \Delta x \frac{\partial \phi_i}{\partial{x}} + \frac{1}{2} \Delta x ^2\frac{\partial^2 \phi_i}{\partial{x^2}} + \mathcal{O}(\Delta x^3)
 \phi_{i-1} = \phi_i - \Delta x \frac{\partial \phi_i}{\partial{x}} + \frac{1}{2} \Delta x ^2\frac{\partial^2 \phi_i}{\partial{x^2}} + \mathcal{O}(\Delta x^3)

Adding the series and solving for \frac{\partial^2 \phi_i}{\partial{x^2}} gives

\frac{\partial^2 \phi_i}{\partial{x^2}}  = \frac{\phi_{i+1} - 2 \phi_i + \phi_{i-1}}{\Delta x^2} + \mathcal{O}(\Delta x^3)

A comprehensive list of these formula, and an interesting discussion is found at [1].

So, having determining the stencil and the level of accuracy we need, we can use Taylor's Theorem to construct suitable finite difference expressions. This is especially important in situations where we might require a special stencil, such as at the edge of the model domain: more on this later. We can also construct finite-difference expressions for grids with non-constant spacing (i.e. where \Delta x varies).

Summary of second-order-accurate FD expressions

First derivatives, second-order accurate, equal spacing


\left.\frac{\partial\phi}{\partial x}\right|_i \approx \frac{\phi_{i+1} - \phi_{i-1}}{2\Delta x}


\left.\frac{\partial\phi}{\partial x}\right|_i \approx \frac{-\phi_{i+2}+4\phi_{i+1}-3\phi_i}{2 \Delta x}
\left.\frac{\partial\phi}{\partial x}\right|_i \approx \frac{3\phi_i-4\phi_{i-1}+\phi_{i-2}}{2 \Delta x}

Second derivatives, second-order accurate, equal spacing


\left.\frac{\partial^2\phi}{\partial x^2}\right|_i \approx \frac{\phi_{i+1}-2\phi_i+\phi_{i-1}}{\Delta x^2}


\left.\frac{\partial^2\phi}{\partial x^2}\right|_i \approx \frac{\phi_{i+2}-2\phi_{i+1}+\phi_{i}}{\Delta x^2}
\left.\frac{\partial^2\phi}{\partial x^2}\right|_i \approx \frac{\phi_{i}-2\phi_{i-1}+\phi_{i-2}}{\Delta x^2}

First derivatives, ?-order accurate, unequal spacing


\left.\frac{\partial\phi}{\partial x}\right|_i \approx \frac{\Delta x_1^2 \phi_{i+1}+(\Delta x_2^2-\Delta x_1^2)\phi_i-\Delta x_2^2 \phi_{i-1}}{\Delta x_1\Delta x_2(\Delta x_1+\Delta x_2)}\qquad\mbox{with}\qquad\Delta x_1 = x_i-x_{i-1} \quad\mbox{and}\quad\Delta x_2 = x_{i+1} - x_i


\left.\frac{\partial\phi}{\partial x}\right|_i \approx \frac{-\Delta x_{1}^2 \phi_{i+2}+(\Delta x_1+\Delta x_2)^2 \phi_{i+1}-(\Delta x_2^2+2\Delta x_1\Delta x_2)\phi_i} {\Delta x_{1}\Delta x_{2}(\Delta x_{1}+\Delta x_{2})}\qquad\mbox{with}\qquad\Delta x_1 = x_{i+1}-x_{i} \quad\mbox{and}\quad\Delta x_2 = x_{i+2} - x_{i+1}
\left.\frac{\partial\phi}{\partial x}\right|_i \approx \frac{(\Delta x_1^2+2\Delta x_1\Delta x_2)\phi_i
-(\Delta x_1+\Delta x_2)^2 \phi_{i-1}
+\Delta x_{2}^2 \phi_{i-2}}
{\Delta x_{1}\Delta x_{2}(\Delta x_{1}+\Delta x_{2})}\qquad\mbox{with}\qquad\Delta x_1 = x_{i-1}-x_{i-2} \quad\mbox{and}\quad\Delta x_2 = x_{i} - x_{i-1}

Second derivatives, ?-order accurate, unequal spacing


\left.\frac{\partial^2\phi}{\partial x^2}\right|_i \approx \frac{2\left[\Delta x_1
    \phi_{i+1}-(\Delta x_1+\Delta x_2)\phi_i +\Delta x_2 \phi_{i-1}\right]}
{\Delta x_{1}\Delta x_{2}(\Delta x_{1}+\Delta x_{2})}\qquad\mbox{with}\qquad\Delta x_1 = x_i-x_{i-1} \quad\mbox{and}\quad\Delta x_2 = x_{i+1} - x_i


\left.\frac{\partial^2\phi}{\partial x^2}\right|_i \approx \frac{2\left[\Delta x_1 \phi_{i+2}
  -(\Delta x_1+\Delta x_2)\phi_{i+1}+\Delta x_2 \phi_i\right]}
{\Delta x_{1}\Delta x_{2}(\Delta x_{1}+\Delta x_{2})}\qquad\mbox{with}\qquad\Delta x_1 = x_{i+1}-x_{i} \quad\mbox{and}\quad\Delta x_2 = x_{i+2} - x_{i+1}
\left.\frac{\partial^2\phi}{\partial x^2}\right|_i \approx \frac{2\left[\Delta x_1 \phi_{i}
  -(\Delta x_1+\Delta x_2)\phi_{i-1}+\Delta x_2 \phi_{i-2}\right]}
{\Delta x_{1}\Delta x_{2}(\Delta x_{1}+\Delta x_{2})}\qquad\mbox{with}\qquad\Delta x_1 = x_{i-1}-x_{i-2} \quad\mbox{and}\quad\Delta x_2 = x_{i} - x_{i-1}


PDEs are typically classified by their characteristics. This scheme is so common that I think it needs to be mentioned. Assuming u_{xy}=u_{yx}, the general second-order PDE in two independent variables has the form

Au_{xx} + Bu_{xt} + Cu_{tt} + \cdots = 0,

where the coefficients A, B, C etc. may depend upon x and t. This form is analogous to the equation for a conic section:

Ax^2 + Bxt + Ct^2 + \cdots = 0.

Just as one classifies conic sections into parabolic, hyperbolic, and elliptic based on the discriminant B^2 - 4AC, the same can be done for a second-order PDE at a given point.

  1. B^2 - 4AC \, < 0 : solutions of elliptic partial differential equation are as smooth as the coefficients allow, within the interior of the region where the equation and solutions are defined. For example, solutions of Laplace's equation are analytic within the domain where they are defined, but solutions may assume boundary values that are not smooth. The motion of a fluid at subsonic speeds can be approximated with elliptic PDEs.
  2. B^2 - 4AC = 0\, : equations that are parabolic partial differential equation at every point can be transformed into a form analogous to the heat equation by a change of independent variables. Solutions smooth out as the transformed time variable increases.
  3. B^2 - 4AC \, > 0  : hyperbolic partial differential equations retain any discontinuities of functions or derivatives in the initial data. An example is the wave equation. The motion of a fluid at supersonic speeds can be approximated with hyperbolic PDEs.


  1. Durran, D. R. (1999) Numerical Methods for Wave Equations in Geophysical Fluid Dynamics. New York: Springer

Other parts of finite differencing

The Lab Classes on Finite Differencing are split up into the following parts: