How Finite-Difference Option-Pricing Algorithms Solve the Black–Scholes PDE: Explicit, Implicit and Crank–Nicolson Grids, Stability, Boundary Conditions and Convergence

Reader question: The Black–Scholes model is usually shown as a closed-form formula for a plain European option. What happens when the payoff, exercise rule or boundary is too complicated for that shortcut?

Finite-difference methods turn an option-pricing partial differential equation into a grid of ordinary arithmetic. Price becomes a set of asset-price nodes, time becomes a set of time steps, derivatives become differences between neighbouring nodes, and the continuous PDE becomes a sequence of linear systems.

This article owns the one-factor PDE discretisation problem: a Black–Scholes-type PDE + terminal payoff + boundary conditions + numerical grid → approximate option values, Greeks, convergence diagnostics and evidence about whether the numerical answer is trustworthy.

It does not own volatility-surface calibration, stochastic-volatility calibration, Monte Carlo path simulation or the derivation of no-arbitrage option pricing itself. Those are separate pages. The purpose here is to show how a continuous pricing equation becomes a computable numerical algorithm.

This is public mathematical and computational education. It is not financial advice, a trading recommendation or a claim that Black–Scholes assumptions describe every market.

1. Start from the PDE

For an option value V(S,t) on an asset price S, constant volatility σ, risk-free rate r and continuous dividend yield q, the Black–Scholes PDE is:

V_t + ½σ²S²V_SS + (r−q)SV_S − rV = 0.

The PDE alone is not enough. We also need a terminal payoff and boundary conditions.

For a European call with strike K and maturity T:

V(S,T) = max(S−K,0).

For a European put:

V(S,T) = max(K−S,0).

The mathematical problem is therefore a backward terminal-value problem: the payoff is known at maturity and we work backward toward today’s value.

2. Reverse time so the algorithm moves forward

It is convenient to define time to maturity:

τ = T − t.

Now τ = 0 at option expiry and increases as the algorithm works backward in calendar time. Since V_t = −V_τ, the PDE can be rewritten as:

V_τ = ½σ²S²V_SS + (r−q)SV_S − rV.

This has the shape of a diffusion–advection–decay equation.

3. Build a grid

Choose a maximum asset price S_max. Divide the interval from 0 to S_max into M steps:

ΔS = S_max/M.

Asset-price nodes are:

S_i = iΔS, i = 0,1,…,M.

Divide time to maturity into N steps:

Δτ = T/N.

Let:

V_iⁿ ≈ V(S_i, nΔτ).

The continuous surface has become a matrix of unknown numbers.

4. Replace derivatives with local differences

At an interior grid point, a central first derivative is:

V_S ≈ (V_i+1 − V_i−1)/(2ΔS).

The central second derivative is:

V_SS ≈ (V_i+1 − 2V_i + V_i−1)/(ΔS)².

These formulas convert the PDE operator at node i into a weighted combination of three neighbouring option values.

That locality creates a tridiagonal matrix, which is one of the main computational advantages of one-dimensional finite-difference pricing.

5. The spatial operator

At node S_i, write the discretised Black–Scholes spatial operator schematically as:

(AV)_i = a_iV_i−1 + b_iV_i + c_iV_i+1.

The coefficients contain:

• diffusion from ½σ²S²V_SS;
• drift from (r−q)SV_S;
• discounting from −rV.

Once A has been assembled, different time-stepping schemes use the same operator in different ways.

6. Explicit finite differences

The explicit Euler scheme uses the known grid at time step n to calculate the next step:

Vⁿ⁺¹ = Vⁿ + Δτ A Vⁿ.

Every new interior value is computed directly from old neighbouring values. No linear solver is required.

The attraction is obvious: the implementation is simple and each step is cheap.

The weakness is stability. If Δτ is too large relative to ΔS, volatility and the location of the grid node, numerical errors can amplify instead of decay. The solution can oscillate or explode even when the mathematical PDE is well behaved.

7. Why explicit stability is not merely a speed setting

For diffusion equations, explicit methods typically require time steps small enough relative to the square of the spatial step. The exact Black–Scholes condition depends on the discretisation and coefficients, but the structural idea is:

finer spatial grid → potentially much smaller stable explicit time step.

Doubling spatial resolution can therefore require far more than doubling the number of time steps.

An explicit grid that “runs” is not necessarily stable. Stability must be tested, not assumed from the absence of a software error.

8. Implicit finite differences

Backward Euler evaluates the spatial operator at the new time level:

Vⁿ⁺¹ = Vⁿ + Δτ A Vⁿ⁺¹.

Rearrange:

(I − Δτ A)Vⁿ⁺¹ = Vⁿ.

The new values are coupled, so each time step requires a tridiagonal linear-system solve.

The benefit is strong stability for the standard parabolic problem. Backward Euler damps high-frequency numerical noise well.

The cost is accuracy: it is first-order in time, so temporal error decreases roughly linearly with Δτ in the asymptotic regime.

9. Crank–Nicolson averages old and new operators

Crank–Nicolson uses the midpoint idea:

(Vⁿ⁺¹ − Vⁿ)/Δτ = ½A Vⁿ + ½A Vⁿ⁺¹.

Rearranging gives:

(I − ½Δτ A)Vⁿ⁺¹ = (I + ½Δτ A)Vⁿ.

For sufficiently smooth solutions, Crank–Nicolson is second-order accurate in time and central differences are second-order in space.

This makes it a common workhorse for one-factor option-pricing PDEs.

But “unconditionally stable” does not mean “immune to ugly numerical behaviour.” A stable scheme can still produce oscillatory local errors near nonsmooth payoffs.

10. One θ-scheme contains all three methods

A useful general form is:

[I − θΔτ A]Vⁿ⁺¹ = [I + (1−θ)Δτ A]Vⁿ.

Then:

• θ = 0: explicit Euler;
• θ = 1: implicit Euler;
• θ = 1/2: Crank–Nicolson.

This shows that the schemes differ mainly in how much of the spatial operator is evaluated at the old versus new time level.

11. Boundary conditions are part of the model

The grid ends at S = 0 and S = S_max, while the theoretical asset-price domain extends to infinity. The PDE cannot solve itself at the truncated boundaries.

For a European call with continuous dividend yield, common asymptotic boundaries are:

V(0,τ) = 0

and for sufficiently large S_max:

V(S_max,τ) ≈ S_maxe^−qτ − Ke^−rτ.

For a European put:

V(0,τ) ≈ Ke^−rτ

and:

V(S_max,τ) ≈ 0.

If S_max is too low, the artificial far boundary contaminates prices closer to the strike.

12. Choosing Smax is a convergence problem

There is no universally correct fixed multiple of the strike. A longer maturity, higher volatility or larger drift uncertainty can push relevant probability mass farther across the grid.

A good implementation increases S_max and checks whether the value at the target spot changes materially.

Falsifier: if moving the far boundary outward changes the reported price beyond tolerance, the original domain was too small.

13. The payoff kink creates a numerical problem

A vanilla option payoff is continuous but not differentiable at S = K. The call payoff changes slope abruptly from 0 to 1. The put changes slope from −1 to 0.

Crank–Nicolson’s high-order accuracy assumes sufficient smoothness. Near maturity, the kink can excite high-frequency grid modes and produce oscillations, especially in Greeks.

The option value may look reasonable while gamma oscillates badly.

14. Rannacher smoothing

A common remedy is to start with a small number of backward-Euler half-steps before switching to Crank–Nicolson.

Backward Euler strongly damps the high-frequency error introduced by the nonsmooth payoff. After the solution has smoothed, Crank–Nicolson recovers higher-order accuracy.

This is often called Rannacher time stepping.

It is a good example of numerical design responding to the specific structure of an option payoff rather than blindly applying one scheme everywhere.

15. The tridiagonal solver

Implicit and Crank–Nicolson time steps require solving:

Bv = d

where B is tridiagonal.

A generic dense linear solve would waste the structure. The Thomas algorithm performs forward elimination and back substitution in:

O(M)

operations for M spatial nodes, rather than the much larger cost of a generic dense solve.

This is why a one-dimensional PDE can be computationally efficient even with thousands of time steps.

16. A minimal algorithm

1. Choose S_max, M and N.
2. Build the asset-price and time grids.
3. Set the terminal payoff at τ = 0.
4. Assemble the finite-difference operator.
5. Apply boundary conditions.
6. For each time step, build the right-hand side.
7. Solve the tridiagonal system if the scheme is implicit or Crank–Nicolson.
8. Apply any exercise constraint for an American option.
9. Move to the next time level.
10. Interpolate the final grid to the actual spot price if it does not lie exactly on a node.
11. Repeat on finer grids to test convergence.

17. American options change the mathematical problem

For an American option, the holder can exercise before expiry. Therefore the option value must satisfy:

V(S,t) ≥ payoff(S).

Where continuing is optimal, the Black–Scholes PDE holds. Where exercise is optimal, the value equals the payoff.

This creates a free-boundary problem or, after discretisation, a linear complementarity problem.

A compact representation is:

V ≥ g,

BV − d ≥ 0,

(V−g)^T(BV−d) = 0,

where g is the exercise payoff at the current time level.

The complementarity condition says that at each node, either the continuation equation is active or the exercise constraint is active.

18. Why “solve then clip once” is not generally enough

A tempting shortcut is:

1. solve the European-style linear system;
2. replace each value by max(value, payoff).

That projection enforces the value floor, but it does not necessarily solve the coupled complementarity system because changing one node can alter the continuation equation at neighbouring nodes.

Methods such as projected SOR, penalty methods and operator-splitting schemes are designed specifically for the American constraint.

19. Non-uniform grids

A uniform grid spends the same resolution far from the strike as near it. But option curvature is often greatest near the strike, barrier or exercise boundary.

A non-uniform transformation can place more nodes in important regions while keeping a large overall domain.

This improves efficiency but complicates derivative formulas because left and right grid spacings differ.

Verification rule: any non-uniform-grid implementation should be benchmarked against the uniform-grid result and an analytic solution where possible.

20. Discrete dividends and barriers

Cash dividends create jumps in the underlying price mapping at known dates. One approach solves the PDE between dividend dates and remaps the option value across the ex-dividend jump.

Barrier options introduce absorbing, rebate or other special boundary conditions at a barrier that may sit inside the normal asset-price domain.

If the barrier lies between grid nodes, naive interpolation can create material error. Grid alignment or specialised barrier corrections may be needed.

This is a general lesson: contract features often appear numerically as boundaries, jumps or constraints.

21. Inputs and outputs

Inputs can include:

• spot price;
• strike;
• maturity;
• interest-rate and dividend curves;
• volatility specification;
• payoff type;
• exercise style;
• barrier or dividend events;
• S_max;
• number of spatial and time steps;
• θ value;
• boundary formulas;
• linear-system solver;
• interpolation rule;
• American-option solver if required.

Outputs can include:

• option value at the requested spot;
• full price grid;
• delta and gamma from grid differences;
• exercise boundary for American options;
• grid-convergence tables;
• boundary-sensitivity diagnostics;
• linear-solver residuals;
• scheme-comparison errors;
• runtime and memory use.

22. Analytic Black–Scholes is the first benchmark

Before using the PDE solver for contracts with no closed form, test it on a European call or put for which the Black–Scholes formula is known.

Compare:

• price error;
• delta error;
• gamma error;
• error across strikes;
• error across maturities;
• error as the grid is refined.

If a solver cannot reproduce the simple analytic case, complexity should not be added.

23. Grid convergence is stronger evidence than one accurate price

A coarse grid can accidentally land close to the analytic answer because different errors cancel.

Run a sequence such as:

(M,N), (2M,2N), (4M,4N).

For a second-order regime, the error should shrink in a pattern consistent with the theoretical order once the grid is sufficiently fine and the solution is smooth enough.

Observed convergence order is itself a diagnostic. If refinement produces erratic error, investigate boundaries, payoff smoothing, interpolation or instability.

24. Shape diagnostics encode no-arbitrage intuition

For standard vanilla calls under ordinary assumptions, the numerical price should be nonnegative, generally increasing in S and convex in S. For puts, monotonicity is reversed while convexity remains.

A grid with negative option values, wild non-monotonicity or negative gamma where convexity should hold is signalling a numerical or modelling problem.

Put–call parity gives another powerful benchmark for European options with consistent conventions.

25. Evidence polarity

Evidence for confidence includes accurate recovery of analytic Black–Scholes prices, smooth convergence under mesh refinement, agreement between implicit and Crank–Nicolson results at fine resolution, small tridiagonal-solver residuals, stable values as S_max moves outward, sensible price monotonicity and convexity, and stable Greeks after appropriate payoff smoothing.

Evidence against confidence includes explicit-scheme explosions, Crank–Nicolson oscillations near payoff kinks, prices sensitive to the far boundary, non-convergent Greeks, American values below intrinsic value, negative prices, violation of known bounds, strong dependence on whether the strike happens to lie exactly on a grid node, or material disagreement with an analytic benchmark.

26. Counterexample: a stable scheme with oscillatory Greeks

Crank–Nicolson can remain mathematically stable while producing alternating grid errors near a nonsmooth payoff. The price curve may look acceptable but gamma can oscillate from positive to negative across adjacent nodes.

Falsifier: inspect gamma and repeat with Rannacher startup steps or a smaller time step. If the oscillation disappears, the issue was numerical smoothness rather than financial structure.

27. Counterexample: a boundary that is too close

A call grid might set S_max = 1.5K for a long-dated high-volatility option. The far-boundary approximation can then influence the value near today’s spot.

Falsifier: double S_max while keeping enough resolution. If the price changes materially, the original truncation was invalid.

28. Counterexample: explicit instability

An explicit solver can produce plausible values for the first few time steps and then erupt into huge positive and negative numbers.

Falsifier: reduce Δτ materially while keeping the spatial grid fixed. If instability disappears, the time step violated the scheme’s stability region.

29. Counterexample: convergence of price but not gamma

Prices integrate local error and can converge faster visually than second derivatives. A trader or risk engine may care more about gamma than about a few cents of price error.

Falsifier: compute Greeks on successively refined grids and under alternative interpolation rules. If gamma is unstable, the solver is not validated for a gamma-sensitive use even if price is acceptable.

30. Counterexample: model error survives perfect numerics

A PDE grid can solve the Black–Scholes equation to machine precision and still produce a poor market price if constant volatility, continuous hedging, frictionless markets or another model assumption is inappropriate.

Falsifier: separate numerical residual from market calibration residual. If the PDE is solved accurately but observed option prices show systematic smile errors, the problem is model specification, not discretisation.

31. Higher dimensions create the curse of dimensionality

With one asset-price dimension and time, the grid is cheap. Add another stochastic factor and the number of nodes multiplies. Three, four or five state variables can make a direct Cartesian grid enormous.

Alternating-direction implicit methods, sparse grids, finite elements and Monte Carlo become attractive as dimensionality grows.

This explains why finite differences are powerful for low-dimensional contracts but not automatically the best choice for every derivative.

32. Alternatives

Closed-form formulas are fastest and most transparent when available.

Binomial and trinomial trees use state lattices and handle early exercise naturally, though convergence can be slower or oscillatory.

Monte Carlo scales better with dimension and path dependence but handles optimal early exercise less directly.

Fourier and FFT methods can be extremely efficient for models with tractable characteristic functions.

Finite elements offer flexible spatial meshes and geometry.

The numerical method should follow the mathematical structure of the contract, not fashion.

33. Connections to the surrounding Bukit Timah Tutor knowledge estate

The analytic baseline is developed in Black–Scholes option-pricing algorithms. This page begins where the closed-form shortcut ends.

Early exercise connects to American-option binomial-tree algorithms, which solve a related stopping problem through a lattice rather than a PDE grid.

High-dimensional pricing connects to Monte Carlo pricing algorithms.

More complex volatility dynamics connect to Heston calibration and to Dupire local-volatility algorithms, both of which can lead to PDE problems with more complicated coefficients or state dimensions.

The full lane is indexed at Finance & Banking Algorithms | Applied Mathematics in Real Financial Systems.

34. What would falsify confidence?

Confidence should be withdrawn if the solver cannot reproduce a known analytic case; if results fail to converge under grid refinement; if the far boundary changes the target price materially; if explicit stability conditions are ignored; if American values violate exercise constraints; if Greeks are unstable; if no-arbitrage shape checks fail; if linear-system residuals are large; or if numerical accuracy is confused with model validity.

35. Verification and update triggers

Preserve the PDE coefficients, grid, boundary rules, θ value, payoff mapping, event dates, solver method, interpolation rule and convergence results for every model version. Unit-test calls, puts, zero-volatility limits where meaningful, near-zero maturity, deep in/out-of-the-money cases and contracts with known alternative benchmarks.

Trigger review when volatility or rate representations change, discrete events are added, a new exercise style is introduced, the grid transformation changes, pricing moves to a new numerical library, Greeks become decision-critical, dimensionality increases, or model calibration creates coefficient regimes not covered by the original stability tests.

36. Primary and high-quality references

• Paolo Brandimarte, Option Pricing by Finite Difference Methods, in Numerical Methods in Finance and Economics, Wiley.
• Daniel J. Duffy, Black–Scholes Finite Differences for the Impatient, treatment of explicit, implicit and Crank–Nicolson schemes.
• Samu Ikonen and Jari Toivanen, Efficient Numerical Methods for Pricing American Options under Stochastic Volatility, for operator-splitting and complementarity methods.
• Robert C. Merton, Theory of Rational Option Pricing, foundational continuous-time option-pricing theory.
• For model validation principles, Federal Reserve Board and OCC, SR 11-7 — Guidance on Model Risk Management.

Educational boundary: Finite differences answer a numerical question: how accurately have we solved the chosen PDE? They do not by themselves answer the separate empirical question: was the chosen financial model appropriate?