Skip to content

How Finite Differences Solve the Black–Scholes PDE

Original question

How can I use finite difference methods to solve option pricing partial differential equations numerically?
Turn time to maturity into an increasing variable, set the payoff at zero time to maturity, approximate the spatial derivatives on a stock grid, and step forward to today's time to maturity. Supply boundary values at both ends of the grid. Use an explicit scheme for transparency or an implicit scheme when its stability properties are better suited to the grid.

Translate the PDE into a time-stepping rule

Write τ=T−t\tau=T-t, use Sj=j ΔSS_j=j\,\Delta S, and assume no dividends. The PDE becomes Vτ=12σ2S2VSS+rSVS−rVV_\tau=\tfrac12\sigma^2S^2V_{SS}+rSV_S-rV. Central differences in spot and a forward difference in τ\tau give the explicit rule
Vjn+1=ajVj−1n+bjVjn+cjVj+1n,V_j^{n+1}=a_jV_{j-1}^n+b_jV_j^n+c_jV_{j+1}^n,

where aj=12Δτ(σ2j2−rj)a_j=\tfrac12\Delta\tau(\sigma^2j^2-rj), cj=12Δτ(σ2j2+rj)c_j=\tfrac12\Delta\tau(\sigma^2j^2+rj), and bj=1−Δτ(σ2j2+r)b_j=1-\Delta\tau(\sigma^2j^2+r). Apply this update only at interior nodes. For a call with zero dividends, use V(τ,0)=0V(\tau,0)=0 and V(τ,Smax⁡)≈Smax⁡−Ke−rτV(\tau,S_{\max})\approx S_{\max}-Ke^{-r\tau}. The second boundary is an asymptotic approximation based on the deep in-the-money call value.

Backward Euler evaluates the spatial operator at the new layer, requiring a tridiagonal solve at each time step. Crank–Nicolson averages the old and new operators. Its higher order on smooth problems can help, although the payoff kink can motivate a few initial backward-Euler substeps to damp oscillations.
The time index nn labels a whole row of prices at one time to maturity. For an explicit step, compute every interior node in the new row from values in the old row. Reusing a value you have already overwritten changes the numerical method. The boundary nodes come from the contract and asymptotic behaviour, so set them at the appropriate time level separately. Refine time while holding the spot grid fixed to diagnose time stepping, then move the upper boundary or refine the spot grid to study spatial error. Changing everything at once obscures the cause.

The two spot derivatives measure different local features of a row of values. Subtracting the two neighbouring values and dividing by 2ΔS2\Delta S estimates the slope. Adding them and subtracting twice the centre, then dividing by (ΔS)2(\Delta S)^2, estimates curvature. Substitute both differences into the PDE before collecting coefficients. The same neighbour can receive a contribution from both the slope term and the curvature term, and their signs need not agree. Checking that decomposition at one interior node is easier than debugging a whole grid of prices.

Verify a grid before trusting the price

Use S0=K=100S_0=K=100, r=5%r=5\%, σ=25%\sigma=25\%, and τ=1\tau=1 year. Take Smax⁡=400S_{\max}=400, 200 spot intervals, and 4,000 time intervals. The stock grid contains 100100 exactly, so there is no interpolation error at the requested price. The Black–Scholes benchmark is about 12.3412.34. Compare the numerical price with that value, halve the spot step and quarter the time step to preserve the leading stability ratio, and then move the boundary outward to identify truncation error separately.

Code the update and the boundary

import numpy as np strike, rate, vol = 100, 0.05, 0.25 m, n, maturity = 200, 4000, 1.0 spots = np.linspace(0, 400, m + 1) dt = maturity / n j = np.arange(1, m) a = 0.5 * dt * (vol**2 * j**2 - rate * j) b = 1 - dt * (vol**2 * j**2 + rate) c = 0.5 * dt * (vol**2 * j**2 + rate * j) values = np.maximum(spots - strike, 0) for step in range(1, n + 1): values[1:-1] = a * values[:-2] + b * values[1:-1] + c * values[2:] values[-1] = spots[-1] - strike * np.exp(-rate * step * dt) print(np.interp(100, spots, values))

What to check

  • Update every interior node from the previous layer. A loop that overwrites the left neighbour during the sweep changes the numerical scheme.
  • Refine time and space together. An explicit solver can become unstable when you refine spot alone.
  • Reapply time-dependent boundaries after each step and benchmark several spot levels.
Use the options pricer to compare the European benchmark with your grid price near the strike.

References and next steps