1. Black-Scholes-Merton PDE to Parabolic Heat Diffusion Transformation

The canonical Black-Scholes-Merton PDE governing option price $V(S, t)$ is:

$$\frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + r S \frac{\partial V}{\partial S} - r V = 0$$

Under logarithmic coordinate transformations $x = \log(S / K)$, dimensionless time $\tau = \frac{1}{2}\sigma^2(T - t)$, and scaling $V(S, t) = K e^{-\alpha x - \beta \tau} u(x, \tau)$ with $\alpha = \frac{1}{2}\left(\frac{2r}{\sigma^2} - 1\right)$ and $\beta = \alpha^2 + \frac{2r}{\sigma^2}$, the PDE transforms into the dimensionless 1D Parabolic Heat Equation:

$$\frac{\partial u}{\partial \tau} = \frac{\partial^2 u}{\partial x^2}$$

2. Crank-Nicolson Discretization for Option Pricing

Discretizing spatial asset price grid $S_j$ and time steps $t_n$, Crank-Nicolson averages spatial central differences implicitly:

$$-\frac{1}{2}\alpha_j V_{j-1}^{n} + (1 - \beta_j) V_j^{n} - \frac{1}{2}\gamma_j V_{j+1}^{n} = \frac{1}{2}\alpha_j V_{j-1}^{n+1} + (1 + \beta_j) V_j^{n+1} + \frac{1}{2}\gamma_j V_{j+1}^{n+1}$$

Achieves unconditional numerical stability with second-order convergence $\mathcal{O}((\Delta S)^2 + (\Delta t)^2)$.

3. Projected Successive Over-Relaxation (PSOR) for American Options

American options present a Free Boundary Linear Complementarity Problem (LCP) due to early exercise intrinsic payoff constraint $V(S, t) \ge h(S)$:

$$\min\left( -\frac{\partial V}{\partial t} - \mathcal{L}_{\text{BS}} V, \; V(S, t) - h(S) \right) = 0$$

PSOR solves the linear system $\mathbf{A}\mathbf{V}^{n} = \mathbf{b}$ iteratively with over-relaxation factor $\omega \in (1, 2)$ while projecting pointwise onto early exercise payoff $h_j$:

$$V_j^{(k+1)} = \max\left( h_j, \; V_j^{(k)} + \omega \frac{b_j - \sum_{m < j} A_{jm} V_m^{(k+1)} - \sum_{m \ge j} A_{jm} V_m^{(k)}}{A_{jj}} \right)$$

4. Python Verification: Crank-Nicolson European Put Option Solver

import numpy as np
from scipy.linalg import solve_banded
from scipy.stats import norm

def verify_crank_nicolson_put(S0: float = 100.0, K: float = 100.0, r: float = 0.05, sigma: float = 0.2, T: float = 1.0):
    """Solves Black-Scholes PDE via Crank-Nicolson finite differences against analytic BS Put."""
    M = 200      # Asset price grid points
    N_t = 200    # Time steps
    S_max = 3.0 * K
    dS = S_max / M
    dt = T / N_t
    
    S = np.linspace(0, S_max, M + 1)
    V = np.maximum(K - S, 0.0)  # Put terminal payoff
    
    j = np.arange(1, M)
    alpha = 0.25 * dt * (sigma**2 * j**2 - r * j)
    beta  = -0.5 * dt * (sigma**2 * j**2 + r)
    gamma = 0.25 * dt * (sigma**2 * j**2 + r * j)
    
    # Tridiagonal matrices M1 * V_new = M2 * V_old
    ab_LHS = np.zeros((3, M - 1))
    ab_LHS[0, 1:] = -gamma[:-1]
    ab_LHS[1, :]  = 1.0 - beta
    ab_LHS[2, :-1] = -alpha[1:]
    
    for _ in range(N_t):
        rhs = alpha * V[:-2] + (1.0 + beta) * V[1:-1] + gamma * V[2:]
        # Boundary conditions at S=0 (V=K*exp(-rt)) and S=S_max (V=0)
        rhs[0] += alpha[0] * 2.0 * K
        V[1:-1] = solve_banded((1, 1), ab_LHS, rhs)
        
    # Interpolate at S0
    pde_price = np.interp(S0, S, V)
    
    # Exact Black-Scholes Put
    d1 = (np.log(S0/K) + (r + 0.5*sigma**2)*T) / (sigma*np.sqrt(T))
    d2 = d1 - sigma*np.sqrt(T)
    bs_put = K * np.exp(-r*T) * norm.cdf(-d2) - S0 * norm.cdf(-d1)
    
    print(f"Exact Black-Scholes Put Price:    {bs_put:.4f}")
    print(f"Crank-Nicolson PDE Solved Price:  {pde_price:.4f}")
    print(f"Absolute Numerical PDE Error:     {abs(pde_price - bs_put):.2e}")

if __name__ == "__main__":
    verify_crank_nicolson_put()