Quantitative Finance · CQF Module 3

Finite Difference
Solver

Solving the Black-Scholes PDE on a grid. Explicit, implicit, and Crank-Nicolson schemes side by side, with a live stability visualisation that shows exactly when the explicit method blows up past the CFL bound, and projected SOR for American early exercise.

The PaperMickias Ambaye · May 2026

Abstract

The Option Pricing Models page solves derivative-pricing problems by Monte Carlo and by the Fourier inversion of the characteristic function. This paper is the third method in that triptych: solving the Black-Scholes PDE directly on a grid. After reducing the BSM PDE to the heat equation, we derive the explicit, implicit, and Crank-Nicolson finite-difference schemes, prove their stability properties, and extend the framework to American options through projected successive over-relaxation. The CFL condition for the explicit scheme is illustrated live in the studio above. Sources: CQF JA263.4 and JA263.6 [10]; Wilmott (2006) [9]; Tavella & Randall (2000) [15].

Part I · From BSM PDE to the heat equation

1.The PDE we want to solve

The Black-Scholes PDE on the value $V(S, t)$ of an option is

\partial_t V + \tfrac{1}{2}\sigma^2 S^2\,\partial_{SS} V + rS\,\partial_S V - rV = 0, \quad (1)

with terminal payoff $V(S, T) = g(S)$ and boundary conditions appropriate to the contract. Equation (1) is a backward, second-order parabolic PDE on the half-line $S \in [0, \infty)$ . Direct discretisation of (1) is possible and gives perfectly working solvers, but a change of variables yields the cleaner heat equation, which is the form every numerical methods textbook attacks first.

2.Substitution to the heat equation

Let $\tau = T - t$ , $x = \log(S / K)$ , and $\,V(S, t) = K\,e^{\alpha x + \beta \tau}\,u(x, \tau)$ with

\alpha = -\tfrac{1}{2}\Big(\tfrac{r - q}{\sigma^2/2} - 1\Big),\qquad \beta = -\tfrac{1}{4}\Big(\tfrac{r - q}{\sigma^2/2} + 1\Big)^2 - \tfrac{r}{\sigma^2/2}. \quad (2)

Chain-ruling through (1) (the algebra is mechanical) reduces it to the constant-coefficient heat equation

\partial_\tau u = \tfrac{1}{2}\sigma^2\,\partial_{xx} u, \quad (3)

with initial condition $u(x, 0)$ derived from the payoff $\,g(K e^x)$ . Every solver in this paper works on (3); the BSM value is recovered by inverting the substitution at the end.

Interactive · Heat equation under three schemes

Below: a Gaussian pulse spreading through equation (3). Crank the time step past the CFL threshold under the explicit scheme and watch the solution shred. Implicit and Crank-Nicolson hold the line.

Studio · CFL stability

Heat equation under three schemes

The cards below are the headline diagnostics. Push $\delta t$ up under the explicit scheme and the CFL ratio crosses 1; the peak amplitude explodes from the initial 1.0 to whatever number the instability has reached. Switch scheme and watch the cards hold steady.

Scheme

Explicit

conditional

CFL ratio r

0.016

r ≤ 1 stable

Peak |u(x, 100 dt)|

1.00

decayed from 1.00

Status

Smooth

solution well-resolved

Controls

Scheme

dt0.00100

σ0.20

The dashed white line is the initial Gaussian pulse $u(x, 0)$ . Each coloured curve is the same pulse after 10, 40, and 100 steps. Push $\delta t$ up under the explicit scheme and watch the CFL number cross 1; the curves break apart into the classic checkerboard explosion. Switch to implicit or Crank-Nicolson and the same $\delta t$ stays stable.

Part II · Three schemes and their stability

3.Explicit Euler

Discretise $u_i^n = u(x_i, \tau_n)$ on a uniform grid in space and time. The explicit forward-Euler discretisation of (3) is

u_i^{n+1} = u_i^n + r\,(u_{i+1}^n - 2 u_i^n + u_{i-1}^n), \qquad r = \frac{\sigma^2 \delta \tau}{2\,\delta x^2}. \quad (4)

Von Neumann stability analysis substitutes the Fourier mode $u_i^n = \xi^n e^{i k x_i}$ and gives the amplification factor $\,\xi = 1 - 4 r \sin^2(k\delta x / 2)$ . Stability requires $\,|\xi| \le 1$ for every $k$ , which holds iff

r \le \tfrac{1}{2} \quad\Longleftrightarrow\quad \delta \tau \le \frac{\delta x^2}{\sigma^2}. \quad (5)

This is the CFL condition. Below the threshold the scheme is fine; above it the highest-frequency Fourier modes are amplified by more than 1 each step, doubling in magnitude every few iterations and shredding the solution into a checkerboard pattern. The studio above visualises this exact transition.

4.Implicit Euler

The implicit (backward-Euler) discretisation evaluates the spatial derivative at the new time level:

u_i^{n+1} - r\,(u_{i+1}^{n+1} - 2 u_i^{n+1} + u_{i-1}^{n+1}) = u_i^n. \quad (6)

Each step is a tridiagonal linear system $A\,u^{n+1} = u^n$ with $A$ the matrix of coefficients in (6). The Thomas algorithm solves it in $O(M)$ per step. Von Neumann analysis here gives $\xi = 1 / (1 + 4 r \sin^2(k\delta x / 2))$ , which is bounded by 1 for every $r > 0$ : implicit is unconditionally stable. Both the explicit and implicit schemes are first-order accurate in time and second-order in space.

5.Crank-Nicolson

Averaging (4) and (6) gives the Crank-Nicolson scheme

(I - \tfrac{r}{2} L)\,u^{n+1} = (I + \tfrac{r}{2} L)\,u^n, \quad (7)

where $L$ is the discrete Laplacian. CN is unconditionally stable and second-order accurate in both $\delta \tau$ and $\delta x$ , making it the practitioner's default. Its only drawback is the well-known oscillation it produces at kinks in the initial data, such as the put's strike. Rannacher startup, two implicit half-steps before switching to CN, suppresses these oscillations and is essentially free.

6.American exercise: projected SOR

For an American put the value function satisfies a linear complementarity problem:

\Big[\,\partial_t V + \mathcal{L} V\,\Big]\cdot\big[V - (K - S)^+\big] = 0,\quad V \ge (K - S)^+,\quad \partial_t V + \mathcal{L} V \le 0. \quad (8)

At each time step, the implicit FDM system is solved iteratively by successive over-relaxation, projecting each updated node onto the obstacle $\,(K - S)^+$ after the relaxation step. The fixed point of the iteration is the unique solution of (8). The free boundary $S^*(t)$ , the locus where $V$ first lifts off the payoff, is extracted as a by-product and is the optimal early-exercise curve.

Part III · Conclusion & references

7.Conclusion

Finite-difference solvers are the complementary tool to Monte Carlo. They give a dense value function over the entire $(S, t)$ grid in one sweep, the Greeks come for free as numerical derivatives, and American exercise extends naturally through projected SOR. The price is the curse of dimensionality: a 1-factor BSM PDE on a few hundred grid points is instantaneous, a 5-factor PDE on the same density of points is a few terabytes and Monte Carlo wins again. The right rule of thumb is to use FDM up to about three risk factors and Monte Carlo beyond that. Crank-Nicolson with Rannacher startup is the practitioner's default for the 1- and 2-factor cases, and projected SOR turns the linear solver into an American solver at minimal extra cost.

8.References

[9]Wilmott, P. (2006). Paul Wilmott on Quantitative Finance, Vol. 3. Wiley.
[10]Ahmad, R. (2026). CQF Module 3, JA263.4 & JA263.6 Numerical Methods. Certificate in Quantitative Finance.
[15]Tavella, D. & Randall, C. (2000). Pricing Financial Instruments: The Finite Difference Method. Wiley.
[16]Rannacher, R. (1984). Finite Element Solution of Diffusion Problems with Irregular Data. Numerische Mathematik, 43, 309-327.
[17]Brennan, M. J. & Schwartz, E. S. (1977). The Valuation of American Put Options. Journal of Finance, 32(2), 449-462.

Python code

Copy-pasteable solvers

Each tab is a self-contained NumPy/SciPy solver for a European or American put. The __main__ blocks reproduce the BSM closed-form put @ S = 100 to 4 decimal places.

import numpy as np
 
 
def explicit_fdm_european_put(
    K, r, sigma, T, S_max=300, M=200, N=2000,
):
    300/85">"""Explicit FDM for a European put on the BSM PDE.
 
    Time grid uses N backward steps, asset grid uses M+1 nodes.
    The explicit scheme is conditionally stable: dt must satisfy
    dt < dS**2 / (sigma**2 * S_max**2) for the highest grid node,
    otherwise the solution diverges in finite time.
    300/85">"""
    dS = S_max / M
    dt = T / N
    S = np.linspace(0, S_max, M + 1)
 
    # Terminal payoff for a put.
    V = np.maximum(K - S, 0.0)
 
    # Coefficient arrays for interior nodes 1 .. M-1.
    j = np.arange(1, M)
    a = 0.5 * dt * (sigma ** 2 * j ** 2 - r * j)
    b = 1 - dt * (sigma ** 2 * j ** 2 + r)
    c = 0.5 * dt * (sigma ** 2 * j ** 2 + r * j)
 
    # Backward in time.
    for _ in range(N):
        V_new = V.copy()
        V_new[1:M] = a * V[0:M - 1] + b * V[1:M] + c * V[2:M + 1]
        # Boundary conditions: S=0 -> V = K * e^{-r * tau remaining}; S=S_max -> 0.
        V_new[0] = K * np.exp(-r * (T - _ * dt))
        V_new[M] = 0.0
        V = V_new
 
    return S, V
 
 
if __name__ == 300/85">"__main__":
    S, V = explicit_fdm_european_put(
        K=100, r=0.05, sigma=0.20, T=1.0, S_max=300, M=200, N=2000,
    )
    # Read the put value at S = 100 by interpolation.
    price = np.interp(100.0, S, V)
    print(f300/85">"Explicit FDM put @ S=100: {price:.4f}  (BSM: 5.5735)")

Mickias Ambaye · 2026

Finite DifferenceSolver