Finite Difference — Explicit Scheme

Motivation: why this matters in quant finance

For European options on a single underlying, Monte Carlo and analytic Black-Scholes are the standard tools. But for early-exercise options (American), path-dependent options (barriers), and exotic local-vol or interest-rate models, the natural numerical method is to solve the Black-Scholes PDE directly on a grid.

Finite difference methods (FDM) discretise both space (asset price) and time (calendar) and march the option value backwards from expiry. They handle early exercise cleanly (compare hold value to exercise value at every node), they handle barriers exactly (kill the value at barrier nodes), and they handle local volatility surfaces transparently.

The simplest FDM is the explicit scheme: at each time step, the value at every grid node is updated by a weighted average of its neighbours from the previous time step. Fast, simple, and shockingly easy to write incorrectly — the explicit scheme has a stability constraint that limits the time step. This note develops the scheme, the stability analysis, and shows why production code uses implicit/Crank-Nicolson instead.

The informal idea

The Black-Scholes PDE for a European option $V(S, t)$ :

\frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + rS\frac{\partial V}{\partial S} - rV = 0,

with terminal condition $V(S, T) =$ payoff.

Discretise: $S_i = i\Delta S$ for $i = 0, 1, \dots, M$ ; $t_n = T - n\Delta t$ for $n = 0, 1, \dots, N$ . Let $V^n_i = V(S_i, t_n)$ .

Approximate derivatives with central differences in $S$ and forward difference in $t$ (looking backward from expiry):

\frac{\partial V}{\partial t} \approx \frac{V^{n+1}_i - V^n_i}{-\Delta t}, \quad \frac{\partial V}{\partial S} \approx \frac{V^n_{i+1} - V^n_{i-1}}{2\Delta S}, \quad \frac{\partial^2 V}{\partial S^2} \approx \frac{V^n_{i+1} - 2V^n_i + V^n_{i-1}}{\Delta S^2}.

Substituting into the PDE and rearranging gives an explicit recurrence:

V^{n+1}_i = a_i V^n_{i-1} + b_i V^n_i + c_i V^n_{i+1},

where $a_i, b_i, c_i$ are coefficients depending on $\sigma, r, i, \Delta t, \Delta S$ . Each grid value at time $n+1$ depends only on three values at time $n$ . Computing one time step is $O(M)$ .

Formal scheme

Defining the index-based form: at grid point $i$ , the coefficients are

\begin{aligned} a_i &= \frac{1}{2}\Delta t\,(\sigma^2 i^2 - r i),\\ b_i &= 1 - \Delta t\,(\sigma^2 i^2 + r),\\ c_i &= \frac{1}{2}\Delta t\,(\sigma^2 i^2 + r i). \end{aligned}

So $V^{n+1}_i = a_i V^n_{i-1} + b_i V^n_i + c_i V^n_{i+1}$ .

Boundary conditions:

$V(S_{\max}, t) \approx S_{\max} - K e^{-r(T-t)}$ for a call (deeply ITM).
$V(0, t) = 0$ for a call.

Terminal condition: $V(S, T) = (S - K)^+$ .

March backwards from $t = T$ to $t = 0$ .

Stability constraint

The explicit scheme is conditionally stable: it diverges unless

\Delta t \le \frac{1}{\sigma^2 M^2}.

Equivalently, $\Delta t / \Delta S^2 \le 1/(\sigma^2 S_{\max}^2)$ (von Neumann analysis).

For a typical setup with $\sigma = 0.2, S_{\max} = 200, M = 200$ : $\Delta t \le 1/(0.04 \cdot 40000) = 6.25 \times 10^{-4}$ . For $T = 1$ , that's $\sim 1600$ time steps just to stay stable. Halving $\Delta S$ requires quadrupling $\Delta t$ — fine grids in $S$ force tiny time steps.

This is the CFL (Courant-Friedrichs-Lewy) condition in disguise.

Algorithm and example

import numpy as np

def fdm_explicit_call(S0, K, T, r, sigma, M=200, S_max=300, theta=1.0):
    """Explicit FDM for European call. Returns V(S0, 0)."""
    dS = S_max / M
    # CFL-safe dt
    dt = 0.9 * 1.0 / (sigma**2 * M**2)
    N = int(np.ceil(T / dt))
    dt = T / N  # adjust for clean stepping
    
    # Grid
    S = np.linspace(0, S_max, M+1)
    V = np.maximum(S - K, 0)  # terminal condition
    
    # Coefficients (vectorised)
    i = np.arange(1, M)  # interior indices
    a = 0.5 * dt * (sigma**2 * i**2 - r * i)
    b = 1 - dt * (sigma**2 * i**2 + r)
    c = 0.5 * dt * (sigma**2 * i**2 + r * i)
    
    for n 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]
        V_new[0] = 0  # boundary at S = 0
        V_new[M] = S_max - K * np.exp(-r * (T - (n+1)*dt))  # boundary at S_max
        V = V_new
    
    return np.interp(S0, S, V)

price = fdm_explicit_call(100, 100, 1.0, 0.05, 0.2)
print(f"FDM explicit price: {price:.4f}")
# FDM explicit price: 10.4530
# Black-Scholes: 10.4506

Within $\$ 0.003 $of the closed form — good for$ M = 200$.

Key properties

Conditionally stable. Need $\Delta t = O(\Delta S^2)$ for stability — a strong constraint that often dominates compute time.
Easy to code. Three lines per time step. No matrix inversions.
Accuracy: $O(\Delta t) + O(\Delta S^2)$ . First-order in time, second-order in space.
Greeks for free. Once the value grid is computed, $\Delta = \partial V/\partial S$ comes from finite differences on the grid; $\Gamma$ from second differences. No additional simulation needed.
Boundary conditions matter. Wrong $S_{\max}$ truncation, or wrong upper boundary value, can introduce large errors. Truncation reflection: place $S_{\max} \approx 4 S_0$ for typical vol.
American options. Modify the update: $V^{n+1}_i = \max(\text{explicit update}, \text{exercise value})$ . This is a simple obstacle PDE ("LCP").

Worked example: convergence rates

Same call as before, vary grid size:

for M in [50, 100, 200, 400, 800]:
    p = fdm_explicit_call(100, 100, 1.0, 0.05, 0.2, M=M, S_max=400)
    print(f"M={M:>4d}: price = {p:.6f}, error = {abs(p - 10.4506):.6f}")

M=  50: price = 10.4624, error = 0.0118
M= 100: price = 10.4537, error = 0.0031
M= 200: price = 10.4513, error = 0.0007
M= 400: price = 10.4508, error = 0.0002
M= 800: price = 10.4506, error = 0.0000

Halving $\Delta S$ reduces error by $\approx 4\times$ — confirms second-order spatial convergence. Time-step is automatically scaled by CFL.

Common confusions and pitfalls

Forgetting CFL. Setting $\Delta t$ too large produces oscillations that grow exponentially. Result: $V \to \pm\infty$ . Always check stability before running.
Wrong upper boundary. $V(S_{\max}, t)$ is approximated; using $V = 0$ instead of $V \approx S - Ke^{-r(T-t)}$ for a call gives nonsense for $S_0$ close to $S_{\max}$ .
Coordinates in $S$ vs $\ln S$ . The PDE has $S^2 \partial^2 V/\partial S^2$ which makes coefficients depend on $S$ . Transforming to $x = \ln S$ gives a constant-coefficient diffusion equation, which is much easier to discretise. Production code uses log-coordinates.
Greeks from coarse grids. Computing delta from finite differences requires a finer grid than computing the price alone — the price is a smooth function, so coarse grids price well; the gradient amplifies grid noise.
American vs European code reuse. American is European with one extra line: V_new = np.maximum(V_new, payoff) at every step. Don't write two separate routines.
Multi-dimensional curse. For a 2D PDE (basket option, stoch vol), the grid has $M^2$ nodes and the explicit scheme requires $\Delta t = O(\Delta S^2 + \Delta V^2)$ — quickly becomes prohibitive. ADI (alternating direction implicit) or splitting methods are required.

Where this goes next

FDM implicit/Crank-Nicolson — unconditionally stable, dominant in production.
American options via projected SOR or PSOR.
ADI for multi-asset PDEs.