M.8.2 Finite-Difference Methods

Basic Prerequisites: 01 · From Zero and Multivariable Calculus.

---

1. Intuition & Practical Objective

A finite-difference method (FDM) replaces the derivatives in a differential equation by differences on a grid, turning a calculus problem into an algebra problem: a system of linear equations you solve (or march) step by step. It is the workhorse for every one-dimensional initial-boundary-value problem - the heat equation, the Black–Scholes PDE, the diffusion of a bond short rate.

The practical objective of this page is the scheme lookup: given a parabolic PDE, which time-stepping scheme to use, what its order is, and - crucially - whether it is stable.

Two facts carry the whole subject:

1. Consistency + stability $⟺$ convergence. This is the Lax equivalence theorem (Duffy Thm 8.1). Taylor-expand the scheme to check consistency; do a von Neumann (Fourier) test to check stability; convergence then comes free.
2. Stability is a joint property of space and time steps. The explicit scheme is correct and useless outside its CFL bound; the implicit and Crank–Nicolson schemes are unconditionally stable but Crank–Nicolson can ring.

The one-sentence essence. “Replace derivatives by divided differences, then prove the replacement is consistent and stable - stability, not the stencil, decides whether your grid produces the solution or a blow-up.”

---

2. Mathematical Ground Truth & Derivations

Divided differences and their truncation error (Duffy eqs. 6.2–6.10). With $D_0,D_{+},D_{-}$ the centred, forward and backward operators:

$$D_{0}f(a)=\frac{f(a+h)-f(a-h)}{2h},D_{+}f(a)=\frac{f(a+h)-f(a)}{h},D_{-}f(a)=\frac{f(a)-f(a-h)}{h},$$ $$D_{0}f(a)=f^{\prime }(a)+\frac{h^2}{6}\cdot \frac{f^{\prime \prime \prime }({\eta }_{+})+f^{\prime \prime \prime }({\eta }_{-})}{2}=f^{\prime }(a)+O(h^2),D_{\pm }f(a)=f^{\prime }(a)\pm \frac{h}{2}{f}^{\prime \prime }(\eta )=f^{\prime }(a)+O(h),$$ $$D_{+}{D}_{-}f(a)=\frac{f(a-h)-2f(a)+f(a+h)}{h^2}=f^{\prime \prime }(a)+\frac{h^2}{4!}[f^{(4)}({\eta }_{+})+f^{(4)}({\eta }_{-})]=f^{\prime \prime }(a)+O(h^2).$$

So: centred first difference is $O(h^2)$, one-sided is $O(h)$, second difference is $O(h^2)$. (Duffy’s printed (6.10) shows $h^4/4!$ - a typographical error; the Taylor expansion forces $h^2$, as flagged in the corpus.)

The $\theta$-method (Duffy eqs. 6.17–6.19, 7.4–7.7). For a semi-discretised system $\dot{U}=\mathcal{L}U$, weight the spatial operator between the old and new time levels by $\theta$:

$$\frac{U^{n+1}-U^n}{k}=(1-\theta )\mathcal{L}{U}^n+\theta \mathcal{L}{U}^{n+1}⟹[I-k\theta \mathcal{L}]U^{n+1}=[I+k(1-\theta )\mathcal{L}]U^n.$$

$\theta$ on new level	Scheme	Order in time	Stability
$0$	Explicit Euler	$O(k)$	conditional: $\lambda =ak/h^2\le \frac{1}{2}$
$\frac{1}{2}$	Crank–Nicolson	$O(k^2)$	unconditional
$1$	Implicit Euler	$O(k)$	unconditional

(Beware the labelling clash between Duffy §7.3.1 and §7.4.1, where the weight on the new level is written differently - carry one convention and do not mix.)

Consistency and order (Duffy Defs. 8.1, 8.3, 8.4). A scheme is consistent if, for the exact solution $v$,

$$v^{n+1}=Q{v}^n+k{G}^n+k{\tau }^n,\Vert {\tau }^n\Vert \to 0\text{as }h,k\to 0;\text{accurate of order }(p,q)⟺\Vert {\tau }^n\Vert =O(h^p)+O(k^q).$$

${\tau }^n$ is the local truncation error. Lax (Thm 8.1): a consistent two-level scheme for a well-posed linear IVP is convergent iff stable.

Stability - von Neumann/Fourier analysis (Duffy §8.3). Substituting $u_j^n={\gamma }^n{e}^{ij\beta h}$ gives the amplification factor $\rho (\xi )$; stability requires $|\rho (\xi )|\le 1$ for all frequencies $\xi$:

$${\rho }_{\text{expl}}(\xi )=1-4\lambda {\sin }^2\frac{\xi }{2}(\lambda =ak/h^2),{\rho }_{\text{impl}}(\xi )=\frac{1}{1+4\lambda {\sin }^2\frac{\xi }{2}},{\rho }_{\text{CN}}(\xi )=\frac{1-2\lambda {\sin }^2\frac{\xi }{2}}{1+2\lambda {\sin }^2\frac{\xi }{2}}.$$

Explicit Euler needs $\lambda \le \frac{1}{2}$; implicit Euler and Crank–Nicolson are unconditionally stable ($|\rho |<1$ for all $\lambda >0$). Note ${\rho }_{\text{CN}}$ goes negative for large $\lambda {\sin }^2(\xi /2)$ - the source of spurious oscillation.

Stability for matrix schemes - Gerschgorin (Duffy Thm 8.2). For $M{u}^{n+1}=Q{u}^n$, the eigenvalues of $A$ lie in ${\bigcup }_i\{z:|z-a_{ii}|\le {\sum }_{j\ne i}|a_{ij}|\}$; Corollary 8.1 gives $\rho (A)\le {\max }_i{\sum }_j|a_{ij}|$. For tridiagonal Toeplitz systems the spectrum is closed-form (Duffy eqs. 7.8, 8.51): ${\lambda }_j=b+2\sqrt{ac}\cos \left(\frac{j\pi }{n+1}\right)$ - a complex eigenvalue signals oscillation.

Richardson extrapolation (Duffy eqs. 6.28–6.36). If $U_k=W+Mk+O(k^2)$, then

$$V_{k/2}\equiv 2U_{k/2}-U_k=W+O(k^2),$$

a second-order method built from two first-order implicit-Euler solves - without Crank–Nicolson’s ringing.

Exponential fitting (Duffy eq. 6.54–6.56, 11.17). For the convection–diffusion model, replacing the diffusion coefficient by the fitting factor

$$\rho =\frac{\mu h}{2}\coth \left(\frac{\mu h}{2\sigma }\right)$$

reproduces the exact solution of the model ODE at the grid points and stays monotone for any $h$ - the standard cure for convection-dominated (near-degenerate) problems. Its limit $\sigma \to 0$ is automatic upwinding.

Maximum principle (Duffy Thms 3.1–3.2, discrete Lemma 11.1). If $\mathcal{L}u\le 0$ and $u\ge 0$ on the boundary, then $u\ge 0$ inside. A scheme that preserves this (a positive-type scheme / M-matrix, $a_{ij}\le 0$ for $i\ne j$) cannot produce negative values - the property every option price must obey.

---

3. Computational Implementation - three schemes, one benchmark, one blow-up

Solve $u_t=u_{xx}$ on $[0,1]$ with $u(0,t)=u(1,t)=0$, $u(x,0)=\sin (\pi x)$; exact solution $u(x,t)=\sin (\pi x)e^{-{\pi }^{2}t}$. The midline value at $t=0.1$ is $e^{-{\pi }^2\cdot 0.1}=0.372708$.

Three lessons in one table: (i) at equal $(M,N)$ Crank–Nicolson is more accurate than either Euler (7.5e-4 vs 1.1e-3 / 2.6e-3); (ii) refining CN by $2\times$ in each direction cuts the error by $4\times$ (7.5e-4 → 1.9e-4), confirming $O(h^2)+O(k^2)$; (iii) pushing $\lambda$ just above $\frac{1}{2}$ ($0.6$) blows the solution to $-272$ - the same code, one parameter outside its stability bound.

---

4. Failure Modes & First-Principles Breakdowns

1. CFL violation (explicit scheme). The explicit operator has eigenvalues $1-4\lambda {\sin }^2(\xi /2)$; once $\lambda >1/2$ the highest-frequency mode has $|\rho |>1$ and every time step amplifies it. The grid that was sign-accurate becomes a numerical bomb - measured above, $-272$ and $1.5\times 10^5$.
2. Crank–Nicolson ringing. CN is unconditionally stable ($|\rho |<1$) but ${\rho }_{\text{CN}}<0$ at high frequencies; a discontinuous initial condition (a payoff kink, a barrier) excites those modes and the solution oscillates around the discontinuity forever. Cure: Rannacher (two implicit-Euler steps, then CN) or Richardson-extrapolated implicit Euler (Duffy Ch 33).
3. Consistency without stability is worthless. A scheme can have a perfectly vanishing truncation error and still diverge; Lax’s theorem only routes stable consistent schemes to convergence. Always test $|\rho (\xi )|\le 1$.
4. Wrong-way differences on convection-dominated problems. Centred differences for large $\mu h/\sigma$ give complex eigenvalues ⇒ oscillation; the cell-Péclet condition $h\le 2\sigma /\mu$ (Duffy eq. 7.14) must hold, or use exponential fitting for uniform convergence independent of $h$.
5. Round-off floor. The second-difference formula halves the available digits; below $h\approx 10^{-4}$ the total error grows (page 01). Never refine blindly.
6. Boundary conditions set the global order. A second-order interior stencil with a first-order boundary treatment is only first-order globally; ghost points or a reduced first-order system at the boundary restore the order (Duffy Ch 1).

---

5. References

• Duffy, Finite Difference Methods in Financial Engineering
• Hull, Options, Futures, and Other Derivatives
• Strikwerda, J. C.: Finite Difference Schemes and Partial Differential Equations

---

6. Connected Graph Bridges

• Base: 01 · From Zero · Multivariable Calculus
• Continue: 03 · Monte Carlo · Index Hub
• Cross-link (pricing application): Pricing · Finite Differences - the BSM/ADI/penalty specialisation of exactly these schemes