The θ-family on the log-transformed PDE

The PDE, in log-spot

The Black–Scholes PDE for the option value V(S,t) is not solved in spot. It is solved in x = ln S, with time-to-maturity τ = T - t:

The PDE is written as an operator on the interior nodes, with two boundary nodes imposed from the known asymptotics of the contract. The interior system has n-1 unknowns; the two ends are Dirichlet conditions set before each step.

The θ-family

One scalar parameter, θ, covers three schemes. The time-stepping writes the solution at the new time as a weighted blend of explicit and implicit evaluations of the spatial operator L:

The implementation does not special-case the three schemes. It computes θ from the enum and runs the same step logic for all three:

enum class FdScheme { Explicit, Implicit, CrankNicolson };

inline double theta_of(FdScheme s) noexcept {
    switch (s) {
        case FdScheme::Explicit:        return 0.0;
        case FdScheme::CrankNicolson:   return 0.5;
        case FdScheme::Implicit:        return 1.0;
    }
    return 1.0;
}

The spatial operator coefficients are the same for every scheme, set once per grid:

const double a  = 0.5 * sig * sig;       // diffusion
const double b  = r - q - 0.5 * sig * sig; // convection
const double alpha  = a/(dx*dx) - b/(2*dx);
const double beta_c = -2*a/(dx*dx) - r;
const double gamma  = a/(dx*dx) + b/(2*dx);

The interior row of the operator is then L Vᵢ = α·Vᵢ₋₁ + β·Vᵢ + γ·Vᵢ₊₁, and a θ-step of size h builds the tridiagonal system from -th·h·α, 1 - th·h·β, -th·h·γ on the lower, diagonal, and upper.

Thomas solver

The tridiagonal system is solved in O(n) by the Thomas algorithm — forward elimination, backward substitution. No pivoting:

double beta = diag[0];
rhs[0] /= beta;
for (size_t i = 1; i < n; ++i) {
    scratch[i] = upper[i-1] / beta;
    beta = diag[i] - lower[i] * scratch[i];
    rhs[i] = (rhs[i] - lower[i] * rhs[i-1]) / beta;
}
for (size_t i = n-1; i-- > 0;)
    rhs[i] -= scratch[i+1] * rhs[i+1];

Rannacher start-up

This is the reason the repository exists. Crank–Nicolson is second-order accurate for smooth data. A vanilla payoff is not smooth — it has a kink at the strike. That kink excites high-frequency modes in the grid. Crank– Nicolson’s amplification factor tends to −1 as frequency rises, so those modes are not damped; they alternate sign. The observable result is oscillation near the strike and a measured order closer to 1 than 2.

The fix is two fully implicit half-steps at the start, before handing over to Crank–Nicolson for the rest of the run. Implicit Euler damps high-frequency modes completely, so once they are gone the rest of the run is second order as advertised. The half-steps are not a refinement — they do not improve the accuracy of a scheme that was already working. They are a correctness fix for a scheme that is otherwise quietly wrong on exactly the payoffs anyone would price with it.

FdConfig cfg;
cfg.scheme = FdScheme::CrankNicolson; // default
cfg.rannacher_steps = 2;               // default — two half-steps

// Disable to observe the degradation directly:
cfg.rannacher_steps = 0;

The start-up applies only to Crank–Nicolson. An implicit scheme already damps the modes and gains nothing; the explicit scheme is a different stability problem entirely. The implementation sets rannacher_steps = 0 for anything that is not Crank–Nicolson.

Stability reporting

The explicit scheme is conditionally stable: it requires a·Δτ / Δx² ≤ ½. Rather than let the caller discover instability from a price that has exploded, the implementation computes and reports the ratio and flags whether it is safe:

out.stability_ratio = a * dt / (dx * dx);
out.stability_ok = (cfg.scheme != FdScheme::Explicit)
                   || (out.stability_ratio <= 0.5);

The test suite checks this directly — an explicitly configured unstable grid is flagged, and a stable one gives a small error. The idea is the same one that runs through the library: a scheme that is wrong should say so, not return a number that looks fine.

Boundary conditions

The two ends of the grid are not solved; they are imposed from the known asymptotics of the contract. Deep out-of-the-money, the option is worthless. Deep in-the-money, a call is the discounted forward and a put is the discounted strike less the discounted spot:

// Call
v_lo = 0.0;
v_hi = s_hi * exp(-q*τ) - K * exp(-r*τ);

// Put
v_lo = K * exp(-r*τ) - s_lo * exp(-q*τ);
v_hi = 0.0;

For an American contract the boundary value is the maximum of the asymptotic value and the payoff — the early-exercise premium can make the boundary worth more than the asymptotic formula. The implementation takes that maximum before each step.

The right-hand side of the interior system receives the boundary contributions of the implicit part of the operator, moved across: -th·h·α·v_lo on the first row and -th·h·γ·v_hi on the last.

American projection

American exercise on a finite difference grid is handled by projecting onto the payoff after each time step — taking the maximum of the computed value and the payoff at every node:

if (opt.exercise == Exercise::American)
    for (size_t i = 0; i <= n; ++i)
        v[i] = std::fmax(v[i], opt.payoff(exp(x[i])));

For a European contract the projection is skipped entirely. The step configuration — explicit, implicit, or Crank–Nicolson — is the same; only the American flag changes whether the projection runs.

Source