Finite difference method
| Bengt Fornberg (2011), Scholarpedia, 6(10):9685. | doi:10.4249/scholarpedia.9685 | revision #121475 [link to/cite this article] |
Finite Differences (FD) approximate derivatives by combining nearby function values using a set of weights. Several different algorithms are available for calculating such weights. Important applications (beyond merely approximating derivatives of given functions) include linear multistep methods (LMM) for solving ordinary differential equations (ODEs) and finite difference methods for solving partial differential equations (PDEs).
Contents |
Introduction
The first derivative is mathematically defined as \[\tag{1} f^{\prime }(x)=\lim\limits_{h\rightarrow 0}\dfrac{f(x+h)-f(x)}h\]
cf. Figure 1. Taylor expansion of \(f(x+h)\) shows that \[\tag{2} \dfrac{f(x+h)-f(x)}h=f^\prime (x)+\dfrac{hf^{\prime \prime}(x)}{2!}+\dfrac{h^2f^{\prime \prime \prime}(x)}{3!}+\ldots \,\,\,=f^\prime (x)+O(h^1)\]
i.e. the approximation \[\tag{3} f^\prime (x)\approx \dfrac{f(x+h)-f(x)}h\]
is accurate to first order. The FD weights at the nodes \(x\) and \(x+h\) are in this case [-1 1]\(/h\ .\) The FD stencil can graphically be illustrated as
| \(\bigcirc\) | \(\leftarrow\) | entry for \(f^\prime \) value \(\tag{4} \{1\}\) | |
| \(\blacksquare\) | \(\blacksquare\) | \(\leftarrow\) | entries for \(f \) value \(\left\{-\frac{1}{h}, \frac{1}{h}\right\}\) |
| \(\uparrow\) | \(\uparrow\) | ||
| \(x\) | \(x + h\) | \(\leftarrow\) | spatial locations |
The open circle indicates a typically unknown derivative value, and the filled squares typically known function values.
FD formulas in 1-D
Algorithms for computing FD weights
Direct solution of a linear system
The most natural way to obtain FD weights is to require that the resulting formula becomes exact for as high degree polynomials as possible. For example, the weights \(w_1, w_2, \ldots, w_n\) to use at locations \( x_1, x_2, \ldots, x_n \) for approximating \(\frac{d}{dx}\) at \(x = x_c\) can be found by solving \[\tag{5} \begin{bmatrix} 1 & 1 & \cdots & 1 \\ x_1 & x_2 & \cdots & x_n \\ \vdots & \vdots & \ddots & \vdots \\ x_1^{n-1} & x_2^{n-1} & \cdots & x_n^{n-1} \end{bmatrix} \begin{bmatrix} w_1 \\ w_2 \\ \vdots \\ w_n \end{bmatrix} = \begin{bmatrix} \frac{d}{dx} 1 \vert_{x = x_c} \\ \frac{d}{dx} x \vert_{x = x_c} \\ \vdots \\ \frac{d}{dx} x^{n-1} \vert_{x = x_c} \end{bmatrix} \]
The successive lines of this system enforce that the weight set leads to the correct result for the functions \(1, x, x^2, \ldots, x^{n-1}\) and thus, by linearity, for all polynomials up through degree \(n-1\ .\) This direct approach is very flexible, but not computationally fast (\(O(n^3)\) operations).
Numerically fast and stable algorithms
The following Matlab code (Fornberg 1998) calculates FD weights on arbitrarily spaced 1-D node sets
function c = weights(z,x,m)
% Calculates FD weights. The parameters are:
% z location where approximations are to be accurate,
% x vector with x-coordinates for grid points,
% m highest derivative that we want to find weights for
% c array size m+1,lentgh(x) containing (as output) in
% successive rows the weights for derivatives 0,1,...,m.
n=length(x); c=zeros(m+1,n); c1=1; c4=x(1)-z; c(1,1)=1;
for i=2:n
mn=min(i,m+1); c2=1; c5=c4; c4=x(i)-z;
for j=1:i-1
c3=x(i)-x(j); c2=c2*c3;
if j==i-1
c(2:mn,i)=c1*((1:mn-1)'.*c(1:mn-1,i-1)-c5*c(2:mn,i-1))/c2;
c(1,i)=-c1*c5*c(1,i-1)/c2;
end
c(2:mn,j)=(c4*c(2:mn,j)-(1:mn-1)'.*c(1:mn-1,j))/c3;
c(1,j)=c4*c(1,j)/c3;
end
c1=c2;
end
For example, the statement weights(0,-2:2,6) returns the output
0 0 1.0000 0 0
0.0833 -0.6667 0 0.6667 -0.0833
-0.0833 1.3333 -2.5000 1.3333 -0.0833
-0.5000 1.0000 0 -1.0000 0.5000
1.0000 -4.0000 6.0000 -4.0000 1.0000
0 0 0 0 0
0 0 0 0 0
This shows the optimal weights to be applied to function values at \( x = -2,-1, 0, 1, 2\ ,\) for approximating the 0th up through the 6th derivative at \(x = 0\ .\) The last two lines are all zero, since there exist no formulas for the 5th and the 6th derivative that extend over only 5 node points.
Differentiation matrices
On non-uniform grids of finite width, the derivative approximation at each node point requires a separate set of weights. If every stencil extends over all the node points, algorithms to calculate them can save operations by utilizing the fact that all the stencils are based on the same node set (Weideman and Reddy 2000).
Exact weights by means of symbolic algebra calculation
In the case of equispaced nodes (with separation \(h\)), a very short and general symbolic algebra algorithm was first given in (Fornberg 1998), where a brief derivation can be found. We can ask what the optimal weights are in a FD formula that relates values of \(f^{(m)}(x)\) at some locations with values of \(f(x)\) at other locations, as illustrated below in Figure 2.
Here, the three numbers \(s\ ,\) \(d\ ,\) and \(n\) completely describe the stencil shape. The optimal FD weights can be calculated in just two lines of symbolic algebra code. In Mathematica 7, these two lines are:
t = PadeApproximant[x^s*(Log[x]/h)^m,{x,1,{n,d}}];
CoefficientList[{Denominator[t],Numerator[t]},x]
The two examples below illustrate typical applications of this algorithm:
Example 1. The choice \(s=0, d=2, n=2, m=2\) describes a stencil of the shape
\(\bigcirc\) \(\bigcirc\) \(\bigcirc\) \(\blacksquare\) \(\blacksquare\) \(\blacksquare\)
- for approximating the second derivative (since \(m = 2\)). The algorithm produces the output
- \[\left\{\left\{ \frac{h^2}{12}, \frac{5h^2}{6}, \frac{h^2}{12}\right\},\left\{ 1, -2, 1 \right\}\right\}\ :\]
corresponding to the implicit 4th order accurate formula for the second derivative:
- \[ \frac{1}{12} f^{\prime \prime}(x-h) +\frac{5}{6} f^{\prime \prime}(x) + \frac{1}{12} f^{\prime \prime}(x+h) \approx \frac{1}{h^2}\left\{ f(x-h) - 2f(x) + f(x+h) \right\} \]
Example 2. The choice \(s=-2, d=2, n=1, m=1\) describes a stencil of the shape
\(\bigcirc\) \(\bigcirc\) \(\bigcirc\) \(\blacksquare\) \(\blacksquare\)
- for approximating the first derivative. The output
- \[\left\{\left\{ \frac{5h}{12}, -\frac{4h}{3}, \frac{23h}{12}\right\},\left\{ -1, 1 \right\}\right\}\ :\]
is readily rearranged as
- \[ f(x+h) = f(x) + \frac{h}{12}\left(23 f^\prime (x) - 16 f^\prime (x-h) + 5 f^\prime (x-2h)\right)\ :\]
which we later on this article (and in somewhat different notation) will encounter as the third order Adams-Bashforth method for solving ODEs.
Errors when applying FD formulas to given functions
The application of FD formulas gives rise to two different types of errors:
Truncation errors: The approximation (3) is first order accurate; its error is \(O(h^1)\ .\) In most applications, higher orders (especially above second order) are preferable. One way to introduce pseudospectral (PS) methods is to consider the limit of increasing orders.
Rounding errors: The entries in the numerator of (3) can in standard arithmetic be evaluated to an accuracy of \(10^{-16}\ ,\) i.e. after the division with \(h\ ,\) the rounding error becomes \(O(10^{-16}/h)\ .\)
Total error: The best accuracy is usually obtained when these different error types match. In the present example, \(O(h)\) matches \(O(10^{-16}/h)\) when \(h \approx 10^{-8}\ ,\) producing a total error around \(10^{-8}\ .\) Higher order FD methods fare better in this type of analysis, but higher derivatives make the situation much worse since, for the \(p\)th derivative, the rounding error typically becomes \(O(10^{-16}/h^p)\ .\) FD formulas are for this reason only rarely used for derivatives beyond the third or fourth.
FD formulas for analytic functions: If a function is known to be analytic, and it can be computed also for complex arguments, Cauchy's integral formula leads to FD approximations that do not use values near the point along the real axis but instead on a circle in the complex plane around the point of approximation. In this case, high accuracy does not require h to be very small, and also high order derivatives (say, 50th, or 100th) become numerically available to full machine precision.
FD formulas in higher-D
On Cartesian lattices
On a \((x,y)\)-grid, any mixed derivative, such as \(\frac{\partial^3}{\partial x \partial y^2}\ ,\) is most easily approximated by a stencil that amounts to approximating in the two directions in sequence, for example first \(\frac{\partial}{\partial x}\) and then, on that result, approximate \(\frac{\partial^2}{\partial y^2}\ .\) Just like in the case of analytic differentiation, the result will not depend on the order in which the partial derivatives were carried out.
In meshfree settings
With scattered nodes, one might try to generalize (5) to require the exact result for multivariate polynomials of increasing degrees. This approach fails because many perfectly reasonable node configurations will result in singular linear systems. A much more robust approach is to replace polynomials with radial basis functions (RBFs), for example of Gaussian type. In 1-D and in a certain limit (flat basis functions), this approach reproduces all previous 1-D results, but it generalizes reliably to scattered nodes in any number of dimensions. Applications of this emerging RBF-FD methodology not only include FD methods on irregularly shaped domains, but also on simpler ones such as the surface of a sphere. Without the need for a surface-bound coordinate system, complications such as pole singularities will not arise.
With the approach just outlined (often denoted RBF-FD in the literature), the stencil for approximating a derivative at a node will include a certain number of nearby nodes. The kd-tree algorithm can very quickly identify sets of neighbors for each node from arbitrarily ordered (or entirely unordered) node lists. In particular, there is no need with RBF-FD implementations to create any unstructured grid that connects adjacent nodes to each other, nor to decompose the space into a collection of small separate elements. As a result, RBF-FD codes can be very simple also in highly irregular 3-D geometries.
Applications of FD formulas
The following are three important applications of FD approximations
- ODE solution by linear multistep methods (LMM)
- PDE solutions; lattice-based
- PDE solutions; in meshfree settings.
The last area is under rapid development, and will be described at a later time. The other two are summarized below.
Other Scholarpedia articles provide surveys of a wide range of ODE (ordinary differential equation) solution techniques, both for initial value and for boundary value problems. In our present context, we limit our short ODE comments to LMM, since any one of these methods - of any order of accuracy - in principle only amounts to a quite immediate application of a single FD formula.
LMM for solving ODEs
All the main types of numerical methods for ODE initial value problems generalize immediately from a single scalar first order ODE to systems of any number of coupled such equation. Furthermore, ODEs of higher order can always be rewritten as such first order systems. Hence, we focus here only on solving
\[\tag{6} y^\prime = F(t, y(t)),\]
with initial condition (IC) \(y(t_0) = y_0\ .\)
Since ODE initial value problems often describe how some function evolves with time, we have changed notation somewhat, letting the independent variable be \(t\) instead of \(x\ ,\) and our unknown function be denoted by \(y(t)\ .\)
The simplest numerical ODE scheme is Forward Euler (FE). With \(k\) denoting the time step, (3) leads to the approximate scheme:
\[\tag{7} y(t+k) = y(t) + k y^\prime(t),\]
When applied to the task of time stepping, \(y^\prime(t)\) is then replaced by \(F(t,y(t))\ ,\) and (7) therefore gives a numerical approximation \(y(t+k)\) in terms of \(y(t)\) and \(F(t,y(t))\ .\) The local error (committed each time step) is \(O(k^2)\) and the global error (accumulated when advancing the ODE over a fixed time, with \(k \rightarrow 0\ ,\) i.e. when using \(O(1/k)\) time steps) becomes \(O(k^1)\) — thus FE is described as a first order method.
The three main numerical ODE solution methods (LMM, Runge-Kutta methods, and Taylor methods) all have FE as their simplest case, but then extend in different directions in order to achieve higher orders of accuracy and/or better stability properties. Of the three approaches, only LMM amount to an immediate application of FD approximations.
A LMM is essentially nothing but calling the independent variable t instead of x, and then replacing the FD approximations shown in (4) by the more general case in ( Figure 2). The main classes of LMM methods (Adams–Bashforth, Adams–Moulton, and Backward Differentiation) arise all as special cases of ( Figure 2), with \(m = 1\) and for different choices for \(s\ ,\) \(d\ ,\) and \(n\ .\) This algorithm provides the weights for all these cases. With \(m = 1\) and accuracy of order \(p \geq 1\ :\)
| Adams–Bashforth | \(s = 1 - p,\) | \(d = p - 1,\) | \(n = 1.\) |
| Adams–Moulton | \(s = 2 - p,\) | \(d = p - 1,\) | \(n = 1.\) |
| Backward Differentiation | \(s = p,\) | \(d = 0,\) | \(n = p.\) |
We have already illustrated the FD stencil shape associated with the AB3 scheme (and written down its coefficients explicitly). Thinking again of the \(t\)-axis as going to the right, the stencil shapes for the fourth order accurate AB4, AM4 and BD4 schemes become as follows:
|
AB4
|
AM4
|
BD4
|
PDE solutions; lattice-based
The idea of using FD for approximating PDEs was first proposed by L.F. Richardson (Richardson 1911). In a revolutionary paper, he notes that "Step-by-step arithmetical methods of solving ordinary differential equations have long been employed..." and he then proceeds with generalizing this concept to 2-D. After a flawless equilibrium calculation of the stresses in a cross-section of the first Aswan Dam over the Nile, he looks to the future: "The extension to three variables is, however, perfectly obvious. One has only to let the third variable be represented by the number of the page of a book of tracing paper." Although our computational hardware is now far more powerful than pencil and paper, the basic concept of FD approximations for PDEs remains the same.
Time dependent PDEs
If a PDE is time dependent, \[ \frac{\partial u}{\partial t} = F(u, u_x, u_y, \ldots, \{\mbox{maybe higher derivatives}\})\]
the most straightforward FD approach is to place a Cartesian grid over the spatial domain. Given an approximation for \(u\) at time \(t\ ,\) \(F\) is approximated at all interior nodes (with boundary information incorporated as needed), and an ODE solver is invoked to advance in time the resulting coupled ODE system (featuring a separate ODE for each node point). This general approach is known as the Method of Lines (MOL). Its main advantages include:
- Easy to reach high accuracy (and thereby high computational efficiency) in both space and time,
- Easy to interchange ODE solvers (such as explicit, implicit, different orders, etc.)
- Very reduced need for user input with regard to time step selection, error control, etc.
Below are two illustrations of this type of approximations for the 1-D heat equation
- PDE\[\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2}\ ,\]
- Initial condition (IC)\[u(x,t) =u_0(x)\ ,\] Boundary condition (BC)\[u(0,t) = f(t), y(1,t) = g(t)\ .\]
Example 3: Simplest FD scheme:
- FD approximation\[\tag{8} \frac{u(x,t+k) - u(x,t)}{k} = \frac{u(x-h,t) - 2 u(x,t) + u(x+h,t)}{h^2}\ :\]
Stencil in \((x,t)\)-plane:
\(\Box\) \(\vert\) \(\blacksquare\) — \(\blacksquare\) — \(\blacksquare\)
- This amounts to a MOL approximation using \(\frac{\partial^2 u}{\partial x^2} \approx \frac{u(x-h,t) - 2u(x,t) + u(x+h,t)}{h^2}\) and then time stepping with FE. The accuracy becomes \(O(k) + O(h^2)\ .\)
Stability analysis (von Neumann analysis; strictly applicable only on periodic or infinite domains): At time level \(t\ ,\) let the solution be a combination of modes \(u(x,t) = \sigma^{t/k}e^{i\omega x}\ ,\) where \(\sigma\) is the factor by which this mode grows in amplitude for each time step. Substitution into the difference scheme gives, after some simplifications, \(\sigma=1-4\frac{k}{h^2}\sin^2\left(\frac{\omega h}{2}\right)\ .\) For the solution to stay bounded for all \(\omega\) as \(k, h \rightarrow 0\ ,\) we need \(\vert\sigma\vert \leq 1\ ,\) i.e. the stability condition for (8) thus becomes \(\frac{k}{h^2} \leq \frac{1}{2}\ .\)
Stability analysis via energy methods may also allow boundary conditions and influence of variable coefficients to be incorporated. Still another very powerful (and particularly easy to apply) approach for stability analysis in MOL contexts is to compare the numerical spectrum of a spatial discretization against the stability domain of an ODE solver (a certain fixed region in the complex plane, specific for each ODE solver).
Example 4: Crank-Nicolson approximation:
- FD approximation\[\frac{u(x,t+k)-u(x,t)}{k} = \frac{1}{2} \left[ \frac{u(x-h,t+k)-2u(x,t+k)+u(x+h,t+k)}{h^2} + \frac{u(x-h,t)-2u(x,t)+u(x+h,t)}{h^2}\right]\ :\]
Stencil in \((x,t)\)-plane:
\(\Box\) — \(\Box\) — \(\Box\) \(\vert\) \(\blacksquare\) — \(\blacksquare\) — \(\blacksquare\)
- Again, this can be viewed (and implemented) as a MOL approximation using \(\frac{\partial^2 u}{\partial x^2} \approx \frac{u(x-h,t)-2u(x,t)+u(x+h,t)}{h^2}\) and time stepping with the second order Adams-Moulton scheme. The accuracy becomes \(O(k^2)+O(h^2)\ ,\) and the scheme is unconditionally stable. With all \(u\)-values unknown on the new time level, the scheme requires for each time step the solution of a tridiagonal linear system.
For a PDE with so much dissipation as the heat equation, small numerical errors often get smoothed out quickly, and it may be less critical to use schemes of very high orders of accuracy. This situation is quite different for purely convective PDEs, for which \(\frac{\partial u}{\partial t} + \frac{\partial u}{\partial x} = 0\) is the simplest model equation. Any small errors that are committed (and which are not damped out with dissipation) will persist for all times, and dispersive errors (numerical approximations that cause high Fourier modes to travel with incorrect velocities) often become a serious issue. A common manifestation is that pulses that should travel unchanged instead leave behind trailing wave trains. The use of very high orders of accuracy can be a good remedy in the case of linear PDEs (even for solutions that are not very smooth, since high orders helps to keep also sharp solution features intact over long times, as discussed in Fornberg (1996), Section 4.2). However, in the presence of strong nonlinearities (such as shocks), entirely different dynamics arise, usually requiring completely different methods, such as essentially non-oscillatory (ENO), weighted ENO (WENO), or flux limiter schemes. These are usually based on finite difference or finite volume type approximations; for details, see specialized articles.
Time independent PDEs
We will here consider only the 2-D Poisson equation
- PDE\[\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = f(x,y)\ ,\]
- BC\[u(x,y)\] specified along a closed boundary in the \((x,y)\)-plane.
Example 5: Centered second order finite differences in \(x\) and \(y\ ;\) use same step \(h\) in both directions.
- The FD scheme with the smallest possible stencil becomes:
- FD approximation\[\frac{u(x-h,y)-2u(x,y)+u(x+h,y)}{h^2}+\frac{u(x,y-h)-2u(x,y)+u(x,y+h)}{h^2} = f(x,y)\ :\]
Stencil in \((x,y)\)-plane:
\(\Box\) \(\vert\) \(\Box\) — \(\Box\) — \(\Box\) \(\vert\) \(\Box\)
- Once the boundary condition is incorporated, the numerical solution will amount to solving a usually large coupled sparse linear systems of equations. Although direct (non-iterative) sparse linear solvers can be very efficient (for ex. as implemented in Matlab and available through the "\" operator when the system matrix is declared "sparse"), iterative solvers are often preferred. A wide range of such choices is available, such as Jacobi, Gauss-Seidel, successive overrelaxation (SOR), alternate direction implicit (ADI), multigrid (MG), etc. For more detailed information on these, see specialized Scholarpedia entries.
Both speed and reliability of elliptic solvers generally benefit if the FD approximation is diagonally dominant (which is the case above; the magnitude of coefficient for \(u(x,y)\) is greater or equal to the sum of the magnitudes of the remaining coefficients). This poses somewhat of a problem for immediately applying higher order FD approximations to the second derivatives in the Poisson equation. Various enhancements are available to get around this, such as Collatz' Mehrstellenverfahren, Richardson extrapolation, and deferred correction.
References
- Fornberg, B., Calculation of weights in finite difference formulas, SIAM Rev. 40:685-691, 1998.
- Richardson, L.F., The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam, Phil. Trans. Royal Soc., London 210:307-357, 1911.
- Weideman, J.A.C. and Reddy, S.C., A Matlab differentiation matrix suite, ACM Trans. Math. Software 26:465-519, 2000. [1]
- Fornberg, B., A practical guide to pseudospectral methods, Cambridge University Press, 1996.
Further reading
- Iserles, A., A first course in the numerical analysis of differential equations, Cambridge Texts in Applied Mathematics, Cambridge University Press, 1996.
- LeVeque, R.J., Finite difference methods for ordinary and partial differential equations: Steady-state and time-dependent problems, SIAM, Philadelphia, 2007.
- Morton, K.W. and Mayers, D.F., Numerical solution of partial differential equations, Cambridge University Press, 1994.
External links
See also
Linear multistep methods, Ordinary differential equations, Partial differential equations, Radial basis functions, Initial value problems, Runge-Kutta methods, Method of Lines
