The Black–Scholes PDE, solved four ways

Any of these methods will land near the analytic price. That is a weak test. What separates a correct implementation from a plausible one is the rate.

A scheme with the drift term mis-signed often prices close too, if you only ever run it at one grid size. So nothing here is checked at one grid size. Every method is refined across a range of resolutions, and the observed order is fitted as the slope of log₂(error) against log₂(resolution). A scheme that happens to be right at one resolution will not hold the right slope across all of them.

The resolutions plotted here are exactly those the C++ test suite asserts on. Fitting an order is sensitive to which points you fit it on, so sampling a different range for the picture would put a second, disagreeing set of numbers on this page.

Crank–Nicolson does not deliver its second order

This is the result the repository exists for. Crank–Nicolson is second-order accurate for smooth data, and a vanilla payoff is not smooth — it has a kink at the strike. That kink excites high-frequency modes, and Crank–Nicolson's amplification factor tends to -1 as frequency rises, so those modes are never damped. They only alternate sign.

time steps (log) — error vs the closed form (log)

The middle line is Crank–Nicolson as textbooks present it. It tracks the first-order guide, alongside implicit Euler, rather than the second-order one — and its error is the worst of the three at every resolution tested.

This is Rannacher start-up, and it is not a refinement. It is a correctness fix for a scheme that is otherwise quietly wrong on exactly the payoffs anyone would price with it. Both variants are kept in the test suite, so the failure reproduces, not only the fix:

FdConfig cfg;
cfg.scheme = FdScheme::CrankNicolson;
cfg.rannacher_steps = 0;   // 2 by default

Each method is documented in detail on its own page — choose one below, or read the full comparison table.

Every method, measured the same way

What is implemented

• Lattices.

  Binomial in three parameterisations — Cox-Ross-Rubinstein, Jarrow-Rudd, Tian — and the Boyle trinomial. American exercise by taking the payoff maximum at every node. The tree refuses to price when the time step drives the risk-neutral probability outside [0, 1], rather than returning an arbitrageable number.

• Finite differences.

  The whole theta-family: explicit, Crank-Nicolson, implicit, with Rannacher start-up. Tridiagonal solves by the Thomas algorithm. The explicit scheme reports its own stability ratio, which must not exceed 1/2 — so an unstable configuration announces itself instead of being discovered from a price that has exploded.

• Finite elements.

  Galerkin with P1 hat functions, both consistent and lumped mass matrices, on the same theta time-stepping. Lumping the mass matrix reproduces the finite-difference scheme to 2 x 10^-6 — an equivalence the test suite demonstrates rather than mentions.

• Monte Carlo.

  Exact terminal sampling, Euler-Maruyama and Milstein, with antithetic variates and a control variate on the terminal spot, measured at 7.5x tighter standard error for the same path count. Every result carries a standard error and a 95% interval.

Three decisions worth the words

Everything is solved in log-spot. In S the PDE coefficients carry S and S^2, so truncation error varies across the grid and the matrix must be rebuilt as the grid changes. Substituting x = ln S makes every coefficient constant: one tridiagonal matrix serves every time step, and the scheme is uniformly accurate.

Monte Carlo refuses American exercise. Forward simulation cannot value an optimal stopping problem — that needs Longstaff-Schwartz regression. It throws rather than returning a number that looks like a price.

American finite differences use explicit projection, not a full linear-complementarity solve, so they are first order in time near the free boundary even under Crank-Nicolson. Stated because it is a real limitation, not omitted because it is inconvenient.

Running it

cmake -B build -DCMAKE_BUILD_TYPE=Release && cmake --build build
./build/convergence_tests      # 28 checks
./build/convergence_demo       # prices one contract every way, with errors

Header-only, no dependency beyond a C++20 compiler. Without CMake, build.sh and build.bat do the same job. Every figure and every point on this page is emitted by tools/convergence_study.cpp into convergence.csv and plotted from that file, so nothing shown here is a number typed in by hand.