Implicit Scheme and Crank-Nicolson for the BS PDE

Medium finite-difference-methods · Module 2 · ~30 min

Cell 1 — Shared utilities: BS analytic, Thomas algorithm, θ-scheme solver

Foundation cell: Black-Scholes closed form for benchmarking, the Thomas algorithm tridiagonal solver (O(M) per call), and the general θ-scheme PDE solver on the log-price grid. The same solver handles explicit (θ=0), fully implicit (θ=1), and Crank-Nicolson (θ=0.5) by switching the theta parameter.

import numpy as np
from scipy.stats import norm

# ── Black-Scholes analytic ────────────────────────────────────────────────
def bs_call(S, K, T, r, sigma):
    if T <= 0:
        return max(S - K, 0.0)
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)

# ── Thomas algorithm: solve tridiagonal (lower, main, upper) · x = rhs ──
def thomas_solve(lower, main, upper, rhs):
    """
    Thomas algorithm for tridiagonal system.
    lower[0] is unused (no sub-diagonal on first row).
    upper[-1] is unused (no super-diagonal on last row).
    Complexity: O(n). Modifies copies — does not mutate inputs.
    """
    n = len(main)
    w = main.copy().astype(float)   # modified diagonal
    d = rhs.copy().astype(float)    # modified RHS

# Forward sweep: eliminate lower diagonal
    for i in range(1, n):
        if abs(w[i - 1]) < 1e-300:
            raise ZeroDivisionError(f"Near-zero pivot at row {i-1}")
        m = lower[i] / w[i - 1]
        w[i] -= m * upper[i - 1]
        d[i] -= m * d[i - 1]

# Back substitution
    x = np.zeros(n)
    x[-1] = d[-1] / w[-1]
    for i in range(n - 2, -1, -1):
        x[i] = (d[i] - upper[i] * x[i + 1]) / w[i]
    return x

# ── θ-scheme PDE solver on log-price grid ────────────────────────────────
def fd_theta_call(S0, K, T, r, sigma, M=100, N=100, x_width=4.0,
                  theta=1.0, rannacher_steps=0, payoff_fn=None):
    """
    θ-scheme for the transformed BS PDE on x = ln(S).

Parameters
    ----------
    theta         : 0 = explicit, 1 = fully implicit, 0.5 = Crank-Nicolson
    rannacher_steps : number of initial fully implicit steps before switching to theta
    payoff_fn     : callable(S_nodes) → payoff array; default = call payoff

Returns
    -------
    price   : float  — interpolated V(0, S0)
    V_final : ndarray shape (M+1,) — option value at t=0 for all grid nodes
    x_nodes : ndarray shape (M+1,)
    """
    D = 0.5 * sigma**2          # diffusion coefficient
    c = r - D                   # drift in log-price

# Spatial grid in log-price
    x_min = np.log(K) - x_width * sigma * np.sqrt(T)
    x_max = np.log(K) + x_width * sigma * np.sqrt(T)
    dx = (x_max - x_min) / M
    x_nodes = np.linspace(x_min, x_max, M + 1)
    S_nodes = np.exp(x_nodes)
    dt = T / N

# Dimensionless coupling coefficients
    r_D = D / dx**2             # diffusion number
    r_c = c / (2.0 * dx)       # advection number
    a_alpha = dt * (r_D - r_c) # lower sub-diagonal weight
    a_gamma = dt * (r_D + r_c) # upper super-diagonal weight
    a_beta  = dt * (2.0 * r_D + r)  # main diagonal contribution (includes discount)

# Pre-build LHS tridiagonal arrays for arbitrary theta (saved for each regime)
    def build_lhs(th):
        lower = np.full(M - 1, -th * a_alpha)
        main  = np.full(M - 1, 1.0 + th * a_beta)
        upper = np.full(M - 1, -th * a_gamma)
        return lower, main, upper

def build_rhs_coeffs(th):
        """Coefficients applied to V^n to form the explicit part of RHS."""
        c_lo = (1.0 - th) * a_alpha   # coeff for V_{i-1}^n
        c_mi = 1.0 - (1.0 - th) * a_beta  # coeff for V_i^n
        c_hi = (1.0 - th) * a_gamma   # coeff for V_{i+1}^n
        return c_lo, c_mi, c_hi

# Terminal condition (tau = 0): payoff
    if payoff_fn is None:
        payoff_fn = lambda S: np.maximum(S - K, 0.0)
    V = payoff_fn(S_nodes).astype(float)

# Time-march backward from tau=0 to tau=T
    for n_step in range(N):
        tau_next = (n_step + 1) * dt
        # Effective theta for this step (Rannacher smoothing)
        th = 1.0 if n_step < rannacher_steps else theta

lo, ma, up = build_lhs(th)
        c_lo, c_mi, c_hi = build_rhs_coeffs(th)

# RHS from explicit part: B·V^n
        rhs = (c_lo * V[:-2] + c_mi * V[1:-1] + c_hi * V[2:]).copy()

# Boundary conditions at tau_next
        V_bc_left  = 0.0                                      # call: S → 0 ⟹ V → 0
        V_bc_right = S_nodes[-1] - K * np.exp(-r * tau_next)  # call: S → ∞

# Add boundary corrections to first and last interior rows
        rhs[0]  += th * a_alpha * V_bc_left
        rhs[-1] += th * a_gamma * V_bc_right

# Solve tridiagonal system A·V^{n+1} = rhs
        V_int = thomas_solve(lo, ma, up, rhs)

# Reconstruct full solution
        V = np.empty(M + 1)
        V[0]    = V_bc_left
        V[1:-1] = V_int
        V[-1]   = V_bc_right

# Interpolate to S0
    x0 = np.log(S0)
    price = float(np.interp(x0, x_nodes, V))
    return price, V, x_nodes

# ── Smoke test ────────────────────────────────────────────────────────────
S0, K, T, r, sigma = 100.0, 100.0, 1.0, 0.05, 0.20
bs_ref = bs_call(S0, K, T, r, sigma)
p_impl, _, _ = fd_theta_call(S0, K, T, r, sigma, M=100, N=100, theta=1.0)
p_cn,   _, _ = fd_theta_call(S0, K, T, r, sigma, M=100, N=100, theta=0.5)

print(f"BS analytic   : {bs_ref:.6f}")
print(f"Implicit FD   : {p_impl:.6f}  (error {abs(p_impl - bs_ref):.2e})")
print(f"Crank-Nicolson: {p_cn:.6f}  (error {abs(p_cn - bs_ref):.2e})")
print()
print("Thomas solve pivot check: if no ZeroDivisionError was raised, pivots are safe.")
print("Setup complete.")

Run cell to see output.

Not run

Cell 2 — Pricing accuracy: implicit vs CN vs BS analytic over a spot ladder

Price a European call across a spot ladder S ∈ {80, 90, 100, 110, 120} using both the implicit and CN schemes, and compare to the BS closed form. Confirms that both schemes produce sub-0.001 errors for M=100, N=100 — well inside acceptable production tolerance.

# Accuracy comparison over a spot ladder
K, T, r, sigma = 100.0, 1.0, 0.05, 0.20
spot_ladder = [80, 90, 100, 110, 120]
M, N = 100, 100

print(f"K={K}, T={T}, r={r}, σ={sigma}, M={M}, N={N}\n")
print(f"{'Spot':>6}  {'BS analytic':>12}  {'Implicit':>10}  {'Err_impl':>10}  {'CN':>10}  {'Err_CN':>10}")
print("-" * 68)

accuracy_pass = True
for S in spot_ladder:
    bs_p   = bs_call(S, K, T, r, sigma)
    p_impl = fd_theta_call(S, K, T, r, sigma, M=M, N=N, theta=1.0)[0]
    p_cn   = fd_theta_call(S, K, T, r, sigma, M=M, N=N, theta=0.5)[0]
    e_impl = abs(p_impl - bs_p)
    e_cn   = abs(p_cn   - bs_p)
    flag = "" if max(e_impl, e_cn) < 0.005 else " <-- LARGE"
    if max(e_impl, e_cn) >= 0.005:
        accuracy_pass = False
    print(f"{S:>6.0f}  {bs_p:>12.6f}  {p_impl:>10.6f}  {e_impl:>10.2e}  {p_cn:>10.6f}  {e_cn:>10.2e}{flag}")

print()
print(f"All errors < 0.005: {'PASS' if accuracy_pass else 'FAIL'}")
print()
print("Observation: CN errors are systematically smaller than implicit errors")
print("for the same grid — a consequence of CN's second-order time accuracy.")

Run cell to see output.

Not run

Cell 3 — Convergence order: log-log error vs Δτ (implicit ≈ 1, CN ≈ 2)

Fix the spatial grid (M=200) and vary N = {10, 20, 40, 80, 160, 320}. Plot |price − BS_analytic| vs Δτ = T/N on a log-log scale. The implicit scheme should show slope ≈ 1 (O(Δτ)); Crank-Nicolson should show slope ≈ 2 (O(Δτ²)).

Expected result: When you halve Δτ, the implicit error halves; the CN error quarters. The log-log slopes will be approximately −1 and −2 respectively.

# Convergence order: error vs Δτ
# Fix M=200 (fine spatial grid), vary N

S0, K, T, r, sigma = 100.0, 100.0, 1.0, 0.05, 0.20
M_FIXED = 200
N_values = [10, 20, 40, 80, 160, 320]

bs_ref = bs_call(S0, K, T, r, sigma)

errors_impl = []
errors_cn   = []
dt_values   = []

for N in N_values:
    dt = T / N
    p_impl = fd_theta_call(S0, K, T, r, sigma, M=M_FIXED, N=N, theta=1.0)[0]
    p_cn   = fd_theta_call(S0, K, T, r, sigma, M=M_FIXED, N=N, theta=0.5)[0]
    errors_impl.append(abs(p_impl - bs_ref))
    errors_cn.append(abs(p_cn   - bs_ref))
    dt_values.append(dt)

dt_arr   = np.array(dt_values)
err_impl = np.array(errors_impl)
err_cn   = np.array(errors_cn)

# Estimate convergence slopes via linear regression on log-log
slope_impl = np.polyfit(np.log(dt_arr), np.log(np.maximum(err_impl, 1e-12)), 1)[0]
slope_cn   = np.polyfit(np.log(dt_arr), np.log(np.maximum(err_cn,   1e-12)), 1)[0]

print(f"M (spatial) = {M_FIXED} (fixed)    BS ref price = {bs_ref:.6f}\n")
print(f"{'N':>5}  {'Δτ':>8}  {'Err implicit':>14}  {'Err CN':>14}")
print("-" * 50)
for i, N in enumerate(N_values):
    print(f"{N:>5}  {dt_values[i]:>8.4f}  {errors_impl[i]:>14.4e}  {errors_cn[i]:>14.4e}")

print()
print(f"Estimated convergence order (log-log slope):")
print(f"  Implicit scheme : {slope_impl:.2f}  (expected ≈ 1.00)")
print(f"  Crank-Nicolson  : {slope_cn:.2f}  (expected ≈ 2.00)")
print()

# Pass/fail checks
convergence_pass_impl = 0.6 < slope_impl < 1.4
convergence_pass_cn   = 1.4 < slope_cn   < 2.6

print(f"Implicit slope in [0.6, 1.4]: {'PASS' if convergence_pass_impl else 'FAIL'}")
print(f"CN slope in [1.4, 2.6]:       {'PASS' if convergence_pass_cn   else 'FAIL'}")
print()
print("Ratio test (halving Δτ):")
for i in range(1, len(N_values)):
    ratio_impl = errors_impl[i-1] / max(errors_impl[i], 1e-12)
    ratio_cn   = errors_cn[i-1]   / max(errors_cn[i],   1e-12)
    print(f"  N {N_values[i-1]:>3}→{N_values[i]:>3}: impl ratio={ratio_impl:.2f} (≈2), CN ratio={ratio_cn:.2f} (≈4)")

Run cell to see output.

Not run

Cell 4 — CN oscillations: digital call with and without Rannacher smoothing

Price a digital call (pays $1 if S_T > K) using pure CN and CN with 2 Rannacher smoothing steps. Print the option value at nodes near K and compute the maximum absolute oscillation in delta across nodes near the kink. Demonstrates that Rannacher smoothing suppresses the oscillations.

Expected result: Without Rannacher, delta near the strike oscillates wildly (alternating positive and large-negative values). With 2 Rannacher steps, the oscillations are suppressed and delta is monotonically decreasing from near-zero to zero as S moves away from K.

# Digital call: payoff = 1 if S > K, else 0
# This discontinuous payoff seeds high-frequency modes that CN propagates without damping.

K_dig, T_dig, r_dig, sigma_dig = 100.0, 0.25, 0.05, 0.20
S0_dig = 100.0

def digital_payoff(S_nodes):
    return (S_nodes > K_dig).astype(float)

M_dig, N_dig = 120, 40  # moderate grid to make oscillations visible

# Run 1: Pure CN (theta=0.5, no Rannacher)
_, V_cn, x_cn = fd_theta_call(
    S0_dig, K_dig, T_dig, r_dig, sigma_dig,
    M=M_dig, N=N_dig, theta=0.5, rannacher_steps=0, payoff_fn=digital_payoff
)

# Run 2: CN with Rannacher smoothing (2 implicit steps first)
_, V_ran, x_ran = fd_theta_call(
    S0_dig, K_dig, T_dig, r_dig, sigma_dig,
    M=M_dig, N=N_dig, theta=0.5, rannacher_steps=2, payoff_fn=digital_payoff
)

S_cn  = np.exp(x_cn)
S_ran = np.exp(x_ran)

# Find nodes near the strike (within 3 nodes either side)
idx_K = int(np.argmin(np.abs(S_cn - K_dig)))
window = slice(max(0, idx_K - 4), min(M_dig, idx_K + 5))

print(f"Digital call: K={K_dig}, T={T_dig}, r={r_dig}, sigma={sigma_dig}")
print(f"Grid: M={M_dig}, N={N_dig}\n")
print(f"Option values V near strike (nodes {idx_K-4} to {idx_K+4}):")
print(f"{'Node':>6}  {'Spot':>8}  {'V (pure CN)':>14}  {'V (Rannacher)':>15}")
print("-" * 52)
S_win  = S_cn[window]
V_win  = V_cn[window]
VR_win = V_ran[window]
for j in range(len(S_win)):
    print(f"{'':>6}  {S_win[j]:>8.3f}  {V_win[j]:>14.6f}  {VR_win[j]:>15.6f}")

# Delta (central difference) near the kink
dx_dig = (x_cn[1] - x_cn[0])
delta_cn  = np.diff(V_cn)  / dx_dig / np.exp(x_cn[:-1])   # dV/dS approx
delta_ran = np.diff(V_ran) / dx_dig / np.exp(x_ran[:-1])

# Maximum abs oscillation = max delta - min delta in the near-strike window
d_win = delta_cn[max(0,idx_K-3):idx_K+3]
r_win = delta_ran[max(0,idx_K-3):idx_K+3]
osc_cn  = float(np.max(d_win) - np.min(d_win))
osc_ran = float(np.max(r_win) - np.min(r_win))

print()
print(f"Delta oscillation near strike (max - min in window):")
print(f"  Pure CN     : {osc_cn:.4f}")
print(f"  Rannacher   : {osc_ran:.4f}")
print()

rannacher_pass = osc_ran < 0.5 * osc_cn
print(f"Rannacher reduces oscillation by >50%: {'PASS' if rannacher_pass else 'FAIL'}")
print()
print("Note: for a smooth payoff (e.g., European call) the oscillations are")
print("negligible and Rannacher smoothing makes no visible difference.")

Run cell to see output.

Not run

Cell 5 — Tridiagonal residual check: verify Thomas solve accuracy

Directly verify the Thomas algorithm by constructing a known tridiagonal system, solving it, and computing the residual $\|A\mathbf{x} - \mathbf{b}\|_\infty$. This confirms that the Thomas implementation is correct and the residual is at machine epsilon.

# Thomas algorithm residual check
# Build a diagonally dominant tridiagonal system of size 200 with known solution.

np.random.seed(42)
n = 200
x_exact = np.random.randn(n)

# Construct a diagonally dominant tridiagonal matrix
lower = -0.4 * np.random.rand(n)   # sub-diagonal (lower[0] unused)
upper = -0.3 * np.random.rand(n)   # super-diagonal (upper[-1] unused)
main  = np.abs(lower) + np.abs(upper) + 1.0  # strict diagonal dominance

# Build RHS: b = A · x_exact
b = np.zeros(n)
b[0]    = main[0] * x_exact[0] + upper[0] * x_exact[1]
b[-1]   = lower[-1] * x_exact[-2] + main[-1] * x_exact[-1]
for i in range(1, n - 1):
    b[i] = lower[i] * x_exact[i-1] + main[i] * x_exact[i] + upper[i] * x_exact[i+1]

# Solve and check residual
x_solved = thomas_solve(lower, main, upper, b)

# Recompute A·x_solved manually to get residual
r = np.zeros(n)
r[0]  = main[0] * x_solved[0] + upper[0] * x_solved[1] - b[0]
r[-1] = lower[-1] * x_solved[-2] + main[-1] * x_solved[-1] - b[-1]
for i in range(1, n - 1):
    r[i] = lower[i]*x_solved[i-1] + main[i]*x_solved[i] + upper[i]*x_solved[i+1] - b[i]

residual_inf = np.max(np.abs(r))
solution_err = np.max(np.abs(x_solved - x_exact))

print(f"Thomas algorithm verification (n={n})")
print(f"  ‖Ax_solved − b‖∞  = {residual_inf:.2e}   (should be ≈ machine ε ≈ 2e-16)")
print(f"  ‖x_solved − x_exact‖∞ = {solution_err:.2e}")
print()

thomas_pass = residual_inf < 1e-10
print(f"Residual < 1e-10: {'PASS' if thomas_pass else 'FAIL'}")
print()
print("Now verify: the pre-factorisation trick for constant-coefficient systems.")
print("Build the LHS tridiagonal of the implicit BS scheme once, then time-march.")

# Demonstrate that pre-factorising the tridiagonal saves work
# for the constant-coefficient case.
S0_demo, K_demo = 100.0, 100.0
T_demo, r_demo, sigma_demo = 1.0, 0.05, 0.20
M_demo, N_demo = 100, 200

D_d = 0.5 * sigma_demo**2
c_d = r_demo - D_d
dx_d = 2 * 4.0 * sigma_demo * np.sqrt(T_demo) / M_demo  # approximate
r_D_d = D_d / dx_d**2
r_c_d = c_d / (2.0 * dx_d)
a_alpha_d = (T_demo / N_demo) * (r_D_d - r_c_d)
a_beta_d  = (T_demo / N_demo) * (2.0 * r_D_d + r_demo)
a_gamma_d = (T_demo / N_demo) * (r_D_d + r_c_d)
diag_dominant = (1.0 + a_beta_d) > (abs(a_alpha_d) + abs(a_gamma_d))
print(f"\nImplicit BS tridiagonal diagonal dominance check: {'PASS' if diag_dominant else 'FAIL'}")
print(f"  main = {1.0 + a_beta_d:.6f} > lower + upper = {abs(a_alpha_d) + abs(a_gamma_d):.6f}")

Run cell to see output.

Not run

Cell 6 — Validation summary

Aggregate all checks from this notebook into a final pass/fail report. Run all cells first.

# Final validation summary — run cells 1–5 first
print("=" * 60)
print("  VALIDATION SUMMARY: Implicit Scheme and Crank-Nicolson")
print("=" * 60)

checks = {
    "Pricing accuracy (all spots, both schemes)"  : accuracy_pass,
    "Implicit convergence order ≈ 1"              : convergence_pass_impl,
    "CN convergence order ≈ 2"                    : convergence_pass_cn,
    "Rannacher reduces digital-call oscillations" : rannacher_pass,
    "Thomas algorithm residual at machine ε"      : thomas_pass,
}

all_pass = True
for label, result in checks.items():
    status = "PASS" if result else "FAIL"
    icon   = "✓" if result else "✗"
    if not result:
        all_pass = False
    print(f"  {icon}  {label:<46} {status}")

print("-" * 60)
print(f"  Overall: {'ALL CHECKS PASSED' if all_pass else 'SOME CHECKS FAILED'}")
print("=" * 60)
print()
print("Key facts confirmed by this notebook:")
print("  • Both implicit and CN are unconditionally stable (no CFL constraint).")
print("  • CN achieves O(Δτ²) vs O(Δτ) for the implicit scheme.")
print("  • CN produces oscillations near discontinuous payoffs;")
print("    Rannacher smoothing (2 implicit steps) suppresses them.")
print("  • Thomas algorithm solves the tridiagonal system to machine precision.")

Run cell to see output.

Not run