Finite Difference Methods

Most PDE cannot be solved by formula, and even when they can (the heat equation, Black–Scholes) the formula is often less useful than a number. Finite difference methods replace derivatives by difference quotients on a grid and the PDE by a large system of algebraic equations. They are the oldest and most transparent numerical methods for PDE, they are what most option pricing software does when no formula exists, and they carry the analytic ideas of Course 6 into the discrete world: the maximum principle becomes a condition on the time step, Fourier analysis becomes von Neumann stability analysis, and well-posedness becomes stability.

The central lesson is the Lax equivalence theorem: a consistent scheme converges if and only if it is stable. Consistency (the scheme approximates the PDE) is easy to check by Taylor expansion and almost never fails. Stability (errors do not grow) is the real constraint, and ignoring it gives nonsense: the explicit scheme for the heat equation, run with a time step $5\%$ above its stability limit, turns an option price of $10.45$ into about $300$, and $20\%$ above, into $-3\times10^{12}$ (Exercise 9.13).

By the end of this chapter you will be able to:

• derive the explicit, implicit and Crank–Nicolson schemes for the heat equation and compute their truncation errors;
• prove stability of the explicit scheme for $r = k/h^2 \leq \frac12$ by the discrete maximum principle, and compute amplification factors by von Neumann analysis;
• prove that a consistent, stable scheme converges with the order of its truncation error;
• explain the CFL condition for transport and wave equations, and why upwinding is needed;
• solve tridiagonal systems in linear time, and set up a Black–Scholes solver in log-price with correct boundary conditions;
• describe spectral methods and the FFT.

9.1Forecasting the weather by hand

9.2Grids and truncation error

Consider the heat equation $u_t = u_{xx}$ on $0 < x < 1$, $t > 0$, with $u(0, t) = u(1, t) = 0$ and $u(x, 0) = g(x)$. Choose a space step $h = \frac1J$ and a time step $k$, grid points $x_j = jh$, $t_n = nk$, and approximations $U^n_j \approx u(x_j, t_n)$. The basic difference quotients are

The ratio

will govern stability: the heat equation's scaling (6A.3 The Heat Equation on ℝⁿ) says one unit of time is worth two units of space, and so does the grid.

The explicit (forward Euler) scheme evaluates $\delta^2$ at the old time level:

The truncation error $\tau^n_j$ is what is left when the exact solution is put into the scheme, written as an approximation of the PDE:

A scheme is consistent if $\tau \to 0$ as $h, k \to 0$, and has order $(p, q)$ if $\tau = O(k^p + h^q)$.

If $r = \frac16$, the leading terms cancel: since $u_{tt} = u_{xxxx}$ for the heat equation, $\frac k2u_{tt} - \frac{h^2}{12}u_{xxxx} = \frac{h^2}{2}\big(r - \frac16\big)u_{xxxx} = 0$, and the scheme becomes fourth order in space, a curiosity known since the 1950s.

9.3Stability by the maximum principle

The proof has a structure worth isolating, because it is general: the error obeys the scheme with the truncation error as a source, and stability (the scheme does not amplify) turns a bound on the source into a bound on the solution, accumulated over $n = t_n/k$ steps. Each step contributes $k\tau$, and there are $t_n/k$ of them, so the global error is of the same order as the local truncation error.

If $r > \frac12$ the coefficient $1 - 2r$ is negative, the averaging property is lost, and the scheme is unstable: the grid mode $(-1)^j$ is multiplied by $1 - 4r < -1$ at every step (Exercise 9.6).

9.4Von Neumann analysis

On a periodic grid (or ignoring boundaries), every scheme with constant coefficients maps the discrete Fourier mode $U_j = e^{ij\theta}$ ($-\pi < \theta \leq \pi$) to a multiple $\xi(\theta)e^{ij\theta}$, where the amplification factor $\xi(\theta)$ is found by substitution. For the explicit scheme,

Since a general grid function is a sum of modes (the discrete Fourier transform), and the scheme acts on each independently, the scheme is stable in $\ell^2$ exactly when $|\xi(\theta)| \leq 1 + Ck$ for all $\theta$: the von Neumann condition (for single-step scalar schemes with constant coefficients it is necessary and sufficient; Morton and Mayers, section 2.7). For (E), $\xi$ ranges over $[1 - 4r, 1]$, and $|\xi| \leq 1$ iff $r \leq \frac12$, the same condition the maximum principle gave. The most dangerous mode is $\theta = \pi$, the sawtooth $(-1)^j$, which no smooth solution contains but rounding errors always do.

Implicit schemes. The backward Euler scheme evaluates $\delta^2$ at the new time level,

with $\xi(\theta) = \frac{1}{1 + 4r\sin^2(\theta/2)} \in (0, 1]$: unconditionally stable, for every $r$. It also satisfies the discrete maximum principle for every $r$ (Exercise 9.7). The price is that each step requires solving a linear system, and the truncation error is still $O(k + h^2)$. The Crank–Nicolson scheme (John Crank and Phyllis Nicolson, 1947) averages the two:

It is centred at $t_{n+1/2}$, hence second order in time, $\tau = O(k^2 + h^2)$, and $|\xi| \leq 1$ for all $r$. But as $r\sin^2\frac\theta2 \to \infty$, $\xi \to -1$: high-frequency modes are not damped but flip sign at every step. For smooth data this is harmless; for data with a kink, like an option payoff, it produces slowly decaying oscillations near the kink, visible in the computed delta and gamma. The standard cure (Rolf Rannacher, 1984) is to take the first few steps with backward Euler, which damps high frequencies strongly, and then switch to Crank–Nicolson.

9.5Convergence: the Lax equivalence theorem

Write a linear one-step scheme as $U^{n+1} = S_kU^n$ (for implicit schemes $S_k$ includes the inverse of the implicit matrix), with $S_k$ acting on grid functions with some norm $\|\cdot\|$. The scheme is stable if $\|S_k^n\| \leq C$ for $nk \leq T$, uniformly as $h, k \to 0$ along the chosen refinement path (for example $r$ fixed).

The converse is a consequence of the uniform boundedness principle (4A.2 Baire Category and Its Consequences): if the scheme converged for every initial value but $\|S_k^n\|$ were unbounded, Banach–Steinhaus would produce initial data for which the computed solutions are unbounded, contradicting convergence. So for linear problems the theory is complete: convergence = consistency + stability, and the work is in proving stability.

9.6Solving the implicit equations: the Thomas algorithm

Implicit schemes in one space dimension require, at each time step, the solution of a tridiagonal system

(with $a_1 = c_{J-1} = 0$ after the boundary values are moved to the right side). Gaussian elimination without pivoting, organised as one forward sweep and one backward sweep, solves it in $O(J)$ operations:

This is the Thomas algorithm (Llewellyn Thomas, 1949). It is stable, with no pivoting needed, when the matrix is diagonally dominant, $|b_j| \geq |a_j| + |c_j|$ (Exercise 9.9), which is the case for the implicit heat and Black–Scholes matrices: their diagonal is $1 + 2r$ (plus $rk$) and off-diagonals $-r$ (plus small drift terms). So an implicit step costs about as much as an explicit one, and the restriction $k \leq \frac{h^2}{2}$ disappears. In two or more space dimensions the systems are no longer tridiagonal, and one uses alternating direction implicit (ADI) splittings, which reduce each step to tridiagonal solves in one direction at a time, or iterative solvers.

9.7Transport, waves and the CFL condition

For the transport equation $u_t + au_x = 0$ ($a > 0$), information travels to the right at speed $a$ (6B.1 First-Order Equations and Characteristics). The explicit scheme with a centred space difference,

has $\xi(\theta) = 1 - i\nu\sin\theta$, with $|\xi|^2 = 1 + \nu^2\sin^2\theta > 1$: it is unconditionally unstable, for every time step. The upwind scheme differences in the direction the information comes from,

a weighted average when $0 \leq \nu \leq 1$, hence stable in the maximum norm, and first order accurate.

The CFL condition is necessary, not sufficient (the centred scheme satisfies it and is still unstable). For the wave equation $u_{tt} = c^2u_{xx}$ the standard leapfrog scheme $\frac{U^{n+1}_j - 2U^n_j + U^{n-1}_j}{k^2} = c^2\delta^2U^n_j$ is stable exactly when $c\frac kh \leq 1$ (Exercise 9.11). Note the different scalings: for hyperbolic problems an explicit scheme needs $k = O(h)$, for parabolic problems $k = O(h^2)$, which is why explicit schemes are natural for waves and implicit ones for diffusion.

9.8Spectral methods and the FFT

For smooth periodic problems, finite differences can be replaced by spectral methods: expand the solution in Fourier modes, $u(x, t) = \sum_k\hat u_k(t)e^{ikx}$, differentiate exactly ($\widehat{u_{xx}} = -k^2\hat u_k$), and truncate to $|k| < N/2$. For the heat equation on the circle this solves the problem exactly mode by mode, $\hat u_k(t) = e^{-k^2t}\hat u_k(0)$ (6A.6 Weak Solutions and Elliptic Regularity). For smooth solutions the truncation error decays faster than any power of $N$, exponentially for analytic solutions, compared with the fixed algebraic rates of finite differences. Nonlinear terms are computed pointwise on the grid and transformed back (pseudospectral methods).

What makes this practical is the fast Fourier transform: the discrete Fourier transform $\hat U_k = \sum_{j<N}U_je^{-2\pi ijk/N}$ of $N$ values costs $O(N^2)$ operations directly but only $O(N\log N)$ if $N$ is a power of $2$, by splitting the sum into even and odd $j$ and recursing (Exercise 9.12). James Cooley and John Tukey published the algorithm in 1965 (Gauss had found it around 1805, in unpublished work on asteroid orbits). In finance, the FFT computes option prices for whole ranges of strikes at once from a model's characteristic function (6D.10 Monte Carlo, Fourier and Finite Element Methods).

9.9History

Richardson's numerical weather prediction appeared in 1922 and Courant, Friedrichs and Lewy's paper in 1928. John Crank and Phyllis Nicolson published their scheme in 1947. Von Neumann's stability analysis was developed at Los Alamos in the 1940s and described in print around 1950, the year of the ENIAC forecast. Llewellyn Thomas's tridiagonal algorithm dates from 1949, the Lax–Richtmyer equivalence theorem from 1956, and the Cooley–Tukey FFT from 1965. Brennan and Schwartz brought finite differences to option pricing in 1977–78.

9.10Exercises

Solution

$U^1_j = \varepsilon(-1)^j\big(1 + r(-1 - 2 - 1)\big) = (1 - 4r)\varepsilon(-1)^j = -3\varepsilon(-1)^j$. After $10$ steps the amplitude is $3^{10}\times10^{-16} \approx 5.9\times10^{-12}$; after $40$ steps, $3^{40}\times10^{-16} \approx 1.2\times10^3$. The exact heat flow damps the mode $e^{i\pi x/h}$ by $e^{-\pi^2t/h^2}$, essentially instantly: the scheme turns the most strongly damped mode into the most strongly amplified.

Solution

Let $|U^{n+1}_j|$ be maximal at $j$, say $U^{n+1}_j = M > 0$. Then $(1 + 2r)M = U^n_j + rU^{n+1}_{j+1} + rU^{n+1}_{j-1} \leq U^n_j + 2rM$, so $M \leq U^n_j \leq \max|U^n|$. The case of a negative minimum is symmetric.

Solution

(a) With $s = 2r\sin^2\frac\theta2 \geq 0$, $\xi = \frac{1 - s}{1 + s} \in (-1, 1]$. (b) Crank–Nicolson: $\frac{0.0363}{0.0090} \approx 4.0$, $\frac{0.0090}{0.0022} \approx 4.0$, so order $\log_24 = 2$. Backward Euler: ratios about $2.9$, $2.6$, $2.4$, approaching $2$ (first order) as the $O(k)$ error dominates. (c) Errors $e_{200} \approx -0.002249$, $e_{400} \approx -0.000562$: $\frac{4e_{400} - e_{200}}{3} \approx \frac{-0.002248 + 0.002249}{3}$, which is below $10^{-6}$: a thousandfold improvement at no extra cost.

Solution

(a) Eliminating $U_{j-1}$ from equation $j$ using the reduced equation $U_{j-1} + c'_{j-1}U_j = d'_{j-1}$ gives $(b_j - a_jc'_{j-1})U_j + c_jU_{j+1} = d_j - a_jd'_{j-1}$; divide by the coefficient of $U_j$. (b) $|c'_1| = |c_1|/|b_1| < 1$; if $|c'_{j-1}| < 1$ then $|b_j - a_jc'_{j-1}| \geq |b_j| - |a_j| > |c_j|$, so $|c'_j| < 1$. (c) About $8J$ arithmetic operations: linear in $J$, compared with $O(J^3)$ for a dense solve.

Solution

$S\partial_S = \partial_x$ and $S^2\partial_{SS} = \partial_{xx} - \partial_x$, and $\partial_t = -\partial_\tau$. Substituting gives the stated equation, with constant coefficients. A uniform $x$-grid is geometric in $S$: fine near small prices, coarse far out, matching the scale of relative price changes; and constant coefficients make the stability analysis of this chapter apply directly.

Solution

The roots have product $1$. With $\beta = 1 - 2\nu^2\sin^2\frac\theta2$, they are complex conjugates of modulus $1$ when $|\beta| \leq 1$, and real with one of modulus $> 1$ when $|\beta| > 1$. Since $\beta \leq 1$ always, the condition is $\beta \geq -1$, i.e. $\nu^2\sin^2\frac\theta2 \leq 1$ for all $\theta$, i.e. $\nu \leq 1$.

Solution

$\hat U_k = \sum_{m<M}U_{2m}\omega^{2mk} + \omega^k\sum_{m<M}U_{2m+1}\omega^{2mk}$, and $\omega^2 = e^{-2\pi i/M}$, so the two sums are length-$M$ DFTs (periodic in $k$ with period $M$). Combining them for all $N$ values of $k$ costs $O(N)$. Unrolling $C(N) = 2C(N/2) + cN$ over $\log_2N$ levels gives $cN\log_2N$.

Solution

$k = 1/N$, so $r_{\text{eff}} = \frac{0.02}{0.0004N} = \frac{50}{N}$: $0.455$, $0.5$, $0.526$, $0.625$. The sawtooth factor $1 - 4r_{\text{eff}}$ is $-0.82$, $-1$, $-1.105$, $-1.5$. For $N = 95$ errors grow by about $10\%$ per step, about $1.105^{95} \approx 1.3\times10^4$ overall from tiny rounding and kink errors; for $N = 80$, by $1.5^{80} \approx 1.3\times10^{14}$. That matches the computed failures ($\sim300$ and $\sim3\times10^{12}$) in order of magnitude.

Solution

(a) $W^{n+1}_P = W_P + \frac{ke^{-2W_P}}{h^2}\big(\sum_{Q\sim P}W_Q - 4W_P\big) = (1 - 4\lambda)W_P + \lambda\sum_{Q\sim P}W_Q$ with $\lambda = \frac{ke^{-2W_P}}{h^2}$; the coefficients are non-negative iff $4\lambda \leq 1$. (b) The diffusion coefficient $e^{-2w}$ is the inverse of the metric's conformal factor: as a region shrinks ($w \to -\infty$), the effective diffusion in coordinates speeds up, the allowed $k$ falls like $e^{2\min w}$, and a fixed time step eventually violates stability before the singular time. Simulations therefore adapt $k$ to the current solution (or rescale, or use implicit steps), just as the analysis of singularities rescales by the curvature (11B.1 Ricci Solitons).