Finite Difference — Implicit and Crank-Nicolson

Motivation: why this matters in quant finance

The explicit scheme is conditionally stable: $\Delta t = O(\Delta S^2)$ is required, which forces tiny time steps for fine grids. For a $1000 \times 1000$ grid (typical for Heston or basket-vol PDEs), explicit needs $\sim 10^7$ time steps — utterly infeasible.

Implicit and Crank-Nicolson schemes are unconditionally stable: any

\Delta t

works (subject to accuracy, not stability). The price: each time step solves a tridiagonal linear system,

O(M)

work per step instead of

O(M)

per step but with

O(M^2)

steps.

Net result: implicit FDM with $\Delta t \sim \Delta S$ runs $O(M)$ times faster than explicit at the same spatial accuracy — typically a 100x to 1000x speedup. Crank-Nicolson is second-order accurate in time as well (vs first-order for fully implicit), making it the workhorse method for production option pricing PDEs.

This note develops both, derives their stability properties, and shows the standard Rannacher trick for handling kinky payoffs.

The informal idea

Recall the explicit update for the BS PDE:

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

(Using "later" index

n+1

for closer to $t = 0$ ; we march backward in time.)

Explicit evaluates the spatial derivatives at time

n

(already known); the update is direct.

Implicit evaluates spatial derivatives at time

n+1

(the unknown); the update requires solving a tridiagonal system:

-a_i V^{n+1}_{i-1} + (1 - b_i + 2)V^{n+1}_i - c_i V^{n+1}_{i+1} = V^n_i,

(coefficients adjusted for the BS PDE; I'll write them properly in the next section). This is a tridiagonal linear system, solvable in $O(M)$ via the Thomas algorithm.

Crank-Nicolson averages the two: half the spatial derivative at

n

, half at

n+1

. Symmetric in time, second-order accurate in

\Delta t

, and (in moderate parameter regimes) stable for any

\Delta t

Formal schemes

Define the spatial operator $\mathcal{L}$ such that $\partial V/\partial t = -\mathcal{L} V$ where $\mathcal{L} V = \frac{1}{2}\sigma^2 S^2 V_{SS} + rSV_S - rV$ (note the sign — going backward in time, $\partial V/\partial t = -\mathcal{L} V$ in our convention). On the discrete grid, $\mathcal{L}$ becomes a tridiagonal matrix with rows

(\mathcal{L} V)_i = \alpha_i V_{i-1} + \beta_i V_i + \gamma_i V_{i+1},

with $\alpha_i = \frac{1}{2}(\sigma^2 i^2 - r i)$ , $\beta_i = -(\sigma^2 i^2 + r)$ , $\gamma_i = \frac{1}{2}(\sigma^2 i^2 + r i)$ .

Implicit (backward Euler):

\frac{V^{n+1} - V^n}{\Delta t} = \mathcal{L} V^{n+1} \quad\Rightarrow\quad (I - \Delta t\, \mathcal{L})V^{n+1} = V^n.

Solve the tridiagonal system at each step.

Crank-Nicolson:

\frac{V^{n+1} - V^n}{\Delta t} = \frac{1}{2}(\mathcal{L} V^{n+1} + \mathcal{L} V^n) \quad\Rightarrow\quad \left(I - \frac{\Delta t}{2}\mathcal{L}\right)V^{n+1} = \left(I + \frac{\Delta t}{2}\mathcal{L}\right)V^n.

One tridiagonal solve per step.

Stability and accuracy

Scheme	Time accuracy	Stability	Cost per step
Explicit	$O(\Delta t)$	$\Delta t \le 1/(\sigma^2 M^2)$	$O(M)$
Implicit	$O(\Delta t)$	unconditional	$O(M)$ tridiagonal solve
Crank-Nicolson	$O(\Delta t^2)$	unconditional (in mild regime)	$O(M)$ tridiagonal solve

So Crank-Nicolson is the Pareto-dominant choice: same per-step cost as implicit, second-order time accuracy. With an $O(\Delta t^2) + O(\Delta S^2)$ overall error, the optimal scaling is $\Delta t \sim \Delta S$ — not $\Delta t \sim \Delta S^2$ as in explicit.

Caveat: Crank-Nicolson can produce spurious oscillations near non-smooth initial data (kinky payoffs), because it doesn't damp high-frequency modes. The fix is the Rannacher start: do the first 2-4 time steps fully implicit, then switch to CN. Fully implicit damps oscillations strongly; once they're gone, CN takes over for accuracy.

Algorithm: Crank-Nicolson with Rannacher start

import numpy as np
from scipy.linalg import solve_banded

def fdm_cn(S0, K, T, r, sigma, M=200, S_max=400, N=200, rannacher_steps=4):
    """Crank-Nicolson FDM for European call with Rannacher initial steps."""
    dS = S_max / M
    dt = T / N
    
    S = np.linspace(0, S_max, M+1)
    V = np.maximum(S - K, 0.0)
    
    i = np.arange(1, M)
    alpha = 0.5 * (sigma**2 * i**2 - r * i)
    beta = -(sigma**2 * i**2 + r)
    gamma = 0.5 * (sigma**2 * i**2 + r * i)
    
    def step_implicit(V, dt):
        # (I - dt*L) V_new = V_old
        # tridiagonal: -dt*alpha*V_{i-1} + (1 - dt*beta)*V_i - dt*gamma*V_{i+1} = V_old
        ab = np.zeros((3, M-1))
        ab[0, 1:] = -dt * gamma[:-1]   # superdiag
        ab[1, :] = 1 - dt * beta       # diag
        ab[2, :-1] = -dt * alpha[1:]   # subdiag
        rhs = V[1:M].copy()
        # apply boundary contributions to RHS
        rhs[0] += dt * alpha[0] * V[0]  # V[0]=0 anyway
        rhs[-1] += dt * gamma[-1] * V[M]  # boundary at S_max
        V_new_interior = solve_banded((1, 1), ab, rhs)
        V_new = V.copy()
        V_new[1:M] = V_new_interior
        return V_new
    
    def step_cn(V, dt):
        # (I - dt/2 L) V_new = (I + dt/2 L) V_old
        d = dt / 2
        ab = np.zeros((3, M-1))
        ab[0, 1:] = -d * gamma[:-1]
        ab[1, :] = 1 - d * beta
        ab[2, :-1] = -d * alpha[1:]
        # RHS = (I + d*L) V_old applied to interior
        Lv = alpha * V[:M-1] + beta * V[1:M] + gamma * V[2:M+1]
        rhs = V[1:M] + d * Lv
        # boundary contributions: half from old, half from new boundary
        rhs[0] += d * alpha[0] * V[0]
        rhs[-1] += d * gamma[-1] * V[M]
        V_new_interior = solve_banded((1, 1), ab, rhs)
        V_new = V.copy()
        V_new[1:M] = V_new_interior
        return V_new
    
    for n in range(N):
        # update boundary at S_max for the new time step
        t_new = (n+1) * dt
        if n < rannacher_steps:
            V = step_implicit(V, dt)
        else:
            V = step_cn(V, dt)
        V[0] = 0
        V[M] = S_max - K * np.exp(-r * (T - t_new))
    
    return np.interp(S0, S, V)

price = fdm_cn(100, 100, 1.0, 0.05, 0.2, M=200, N=200)
print(f"CN price: {price:.4f}")  # ~ 10.4506

For $M = 200, N = 200$ : error $< \$ 0.001 $vs Black-Scholes. Compare with explicit at$ M = 200 $which required$ N \sim 1600$ to stay stable — CN is 8x fewer time steps with comparable accuracy.

Key properties

Unconditional stability. Allows $\Delta t \gg \Delta S^2$ , dramatically faster for fine grids.
Crank-Nicolson is second-order in time. Error $O(\Delta t^2 + \Delta S^2)$ , vs $O(\Delta t + \Delta S^2)$ for explicit/implicit.
Tridiagonal solves are fast. Thomas algorithm is $O(M)$ per step. No matrix factorisation, no fill-in.
Smoothing for kinks. Rannacher trick: 2-4 fully implicit steps first to damp the discontinuity at $S = K$ . Without it, CN oscillates near the kink.
American options. Replace each tridiagonal solve with a projected version (PSOR or LCP solver) to enforce the early-exercise constraint at every step.
Multi-dimensional generalisation. ADI (alternating direction implicit) splits 2D into a sequence of 1D tridiagonal solves per direction. Standard for stochastic-vol pricing.

Worked example: convergence comparison

Same call, compare explicit vs implicit vs CN at fixed compute budget:

# Explicit needs N ~ 0.9/(sigma^2 M^2) * T
# CN can use N ~ M

for M in [50, 100, 200]:
    # Explicit
    p_ex = fdm_explicit_call(100, 100, 1, 0.05, 0.2, M, S_max=400)
    # CN with N = M
    p_cn = fdm_cn(100, 100, 1, 0.05, 0.2, M, S_max=400, N=M)
    print(f"M={M:>3d}: explicit={p_ex:.5f}, cn={p_cn:.5f}, BS=10.4506")

M= 50: explicit=10.45824, cn=10.45049, BS=10.4506
M=100: explicit=10.45350, cn=10.45069, BS=10.4506
M=200: explicit=10.45154, cn=10.45062, BS=10.4506

CN is 5-10x more accurate at the same $M$ , while using far fewer time steps.

Common confusions and pitfalls

CN oscillations. Without Rannacher, kinky payoffs (especially digitals) produce spurious oscillations near the strike. Looks like noise on the value surface; can be diagnosed by comparing CN to fully implicit (which doesn't oscillate but is less accurate).
Boundary handling. Tridiagonal systems with non-zero boundary values require adding boundary contributions to the RHS. Easy to forget; easy to get wrong by a factor of 2 in the CN average.
Tridiagonal vs full LU. Always use tridiagonal solvers (scipy.linalg.solve_banded or hand-coded Thomas). Full LU for an $M \times M$ system is $O(M^3)$ — kills performance for any reasonable $M$ .
For local vol or non-uniform grids. The coefficients $\alpha_i, \beta_i, \gamma_i$ depend on $S_i$ — this is fine for tridiagonal solver, but if you're using non-uniform grids, the FD formulas are different. Don't blindly copy uniform-grid formulas.
Implicit can mask kinks. Fully implicit is so heavily damping that it can suppress real features near the kink. Use Rannacher to limit this to 2-4 initial steps, then CN takes over.
Stability "in mild regime." CN is provably unconditionally stable for the heat equation. For the BS PDE with convection ( $rS$ term), it's stable for moderate $\Delta t$ but can become marginally stable for very large $\Delta t$ in convection-dominated regimes (e.g., extremely high $r$ or very long $T$ ). In practice almost always fine.

Where this goes next

American options — projected SOR for the obstacle problem.
ADI for multi-asset PDEs.
Local vol calibration via Dupire — uses the FDM grid directly.