How ADI Algorithms Solve the Two-Dimensional Heston PDE: Mixed Derivatives, Douglas/Craig–Sneyd Splitting, Degenerate Boundaries and Convergence

Reader question: The Heston model has two interacting state variables—asset price and variance—and a mixed derivative created by correlation. How can a numerical algorithm solve that two-dimensional pricing PDE without paying the cost of one enormous implicit matrix solve at every time step?

Alternating Direction Implicit, or ADI, algorithms exploit structure. They split the Heston spatial operator into pieces associated with the asset-price direction, the variance direction and the correlation-driven mixed derivative. Instead of solving one large two-dimensional implicit system, they perform a sequence of much cheaper one-dimensional implicit solves.

This article owns the two-dimensional Heston PDE time-stepping problem: Heston risk-neutral parameters + payoff + spatial grids + boundary rules → a stable numerical pricing surface, convergence diagnostics and evidence about whether the PDE solution is trustworthy.

It does not own Heston parameter calibration, one-dimensional Black–Scholes finite differences, Monte Carlo variance simulation or volatility-surface construction. Those have separate roles in the Bukit Timah Tutor finance-and-banking-algorithms estate.

This is public mathematical and computational education. It is not financial advice, a trading recommendation or a claim that the Heston model is empirically correct for every market.

1. Start from the Heston state equations

Under a risk-neutral measure, a common Heston specification is:

dS = (r−q)S dt + √v S dW_S

dv = κ(θ−v)dt + σ_v√v dW_v

with:

dW_SdW_v = ρ dt.

Here S is the asset price, v is instantaneous variance, κ is mean reversion, θ is long-run variance, σ_v is volatility of variance, ρ is spot–variance correlation, r is the domestic risk-free rate and q is a continuous yield or foreign-rate analogue where appropriate.

The numerical problem begins after these risk-neutral parameters have been specified or calibrated. An exact transform formula may exist for some European claims, but many exercise features, barriers and state-dependent payoffs still motivate PDE solvers.

2. Correlation creates the mixed derivative

For an option value V(S,v,t), the Heston pricing PDE can be written:

V_t + ½vS²V_SS + ρσ_vvS V_Sv + ½σ_v²vV_vv + (r−q)S V_S + κ(θ−v)V_v − rV = 0.

The term:

ρσ_vvS V_Sv

is the key complication. It couples movement in the price direction to movement in the variance direction. If ρ = 0, the cross derivative disappears. When |ρ| is large, the interaction becomes numerically important.

3. A full implicit solve is mathematically simple but computationally expensive

Suppose the grid has M_S price nodes and M_v variance nodes. There are roughly:

M_SM_v

unknown option values per time level.

A fully implicit discretisation couples many of those unknowns. The resulting sparse matrix is far larger and structurally more difficult than the tridiagonal matrix of a one-dimensional Black–Scholes grid.

ADI changes the computational geometry. Rather than invert the full two-dimensional operator, it alternates between directions so that each implicit stage has essentially one-dimensional structure.

4. Split the operator into three parts

After spatial discretisation, write the semi-discrete system schematically as:

U′(τ) = A U(τ),

where U is the vector of grid values and τ is time to maturity.

Split:

A = A₀ + A₁ + A₂.

• A₀ contains the mixed S–v derivative.
• A₁ contains derivatives mainly in the asset-price direction.
• A₂ contains derivatives mainly in the variance direction.

Reaction terms such as −rV can be distributed between directional operators according to the chosen implementation.

The purpose of the split is not merely notation. It determines which pieces can be solved implicitly as narrow-banded one-dimensional systems.

5. Why the mixed term is usually treated explicitly

If the cross derivative were made fully implicit together with both directional terms, the easy tridiagonal structure would be lost.

ADI schemes therefore commonly treat A₀ explicitly while correcting for it in later stages. This preserves cheap directional solves while still allowing second-order temporal accuracy in schemes such as Craig–Sneyd, Modified Craig–Sneyd and Hundsdorfer–Verwer under appropriate conditions.

This trade-off is the centre of the method:

keep expensive coupling explicit; keep stiff one-direction diffusion implicit.

6. Douglas splitting: the simplest ADI idea

A Douglas-style step first makes a predictor using the full operator and then corrects sequentially in each implicit direction.

In schematic form:

Y₀ = Uⁿ + Δτ A Uⁿ

followed by solves resembling:

(I−θΔτA₁)Y₁ = Y₀ − θΔτA₁Uⁿ

(I−θΔτA₂)Uⁿ⁺¹ = Y₁ − θΔτA₂Uⁿ.

The exact implementation depends on the operator convention. The important point is that each solve acts mainly along one coordinate direction.

When a mixed derivative is present, the basic Douglas scheme does not generally deliver the same second-order temporal accuracy as the more elaborate correction schemes. It remains useful as a conceptual baseline and sometimes as a damping stage.

7. Craig–Sneyd adds a mixed-term correction

The Craig–Sneyd scheme performs directional implicit corrections and then adds a correction related to the mixed derivative before another set of directional solves.

The extra stage repairs part of the splitting error caused by treating the mixed term explicitly.

For the classical Craig–Sneyd choice θ = 1/2, the scheme is second-order in time under its standard smoothness assumptions and has well-studied stability properties for two-dimensional convection–diffusion equations with mixed derivatives.

8. Modified Craig–Sneyd generalises the correction

The Modified Craig–Sneyd, or MCS, scheme adjusts both the mixed-term correction and the treatment of the full operator. It is designed to retain second-order temporal accuracy while allowing a wider choice of θ.

In Heston numerical studies, θ = 1/3 is a widely analysed and effective setting. The theoretical stability literature gives important unconditional-stability results for diffusion equations with mixed derivatives at appropriate θ values, while practical Heston implementations still need empirical grid and parameter tests because convection, degeneracy and payoff nonsmoothness matter.

“Unconditionally stable” should never be translated into “any grid is accurate.” Stability prevents error explosion; it does not guarantee small discretisation error.

9. Hundsdorfer–Verwer adds another second-order route

The Hundsdorfer–Verwer, or HV, scheme also uses a predictor–corrector structure with directional implicit solves. It is second-order in time for its intended formulation and has strong stability properties.

In comparative Heston studies, both MCS and HV can perform well across challenging parameter sets. The better choice depends on error tolerance, grid design, correlation, Feller-condition behaviour, payoff type and implementation details.

A robust engine should not choose a scheme by reputation alone. It should benchmark the schemes against trusted reference prices and convergence tests.

10. The spatial grid matters as much as the time scheme

Uniform grids are easy to derive but inefficient. Option curvature is concentrated near economically important regions such as:

• the strike;
• barriers;
• low variance;
• early-exercise boundaries.

A non-uniform grid can concentrate nodes where derivatives change quickly and use fewer nodes elsewhere.

This reduces computational cost for a given accuracy, but finite-difference weights become asymmetric and must be derived for unequal spacing. Any grid transformation therefore needs independent unit tests.

11. The boundary v = 0 is mathematically degenerate

At v = 0, the variance-diffusion coefficient:

½σ_v²v

vanishes. The price-diffusion and mixed-derivative terms involving v also vanish.

The PDE becomes degenerate rather than behaving like an ordinary interior two-dimensional diffusion equation. A common numerical treatment is to apply the reduced PDE itself at the boundary using one-sided differences where needed, rather than impose an arbitrary extra boundary value.

This is especially important when the Feller condition:

2κθ ≥ σ_v²

is violated. Violation does not make the Heston model unusable, but it makes the near-zero-variance region more numerically consequential.

12. The Feller condition is not a simple pass/fail test for the PDE solver

The Feller condition is sufficient for the CIR-type variance process to stay strictly positive rather than hit zero. Market calibrations can violate it.

A good solver therefore needs to work in both regimes:

• Feller satisfied;
• Feller violated.

Testing only convenient parameter sets creates false confidence. Published Heston ADI studies explicitly include both cases because numerical difficulty can increase when variance spends more time near zero.

13. The far variance boundary is artificial

The theoretical variance state extends without a fixed upper limit, but a computer grid must truncate it at v_max.

Possible numerical boundary choices include zero-gradient assumptions or asymptotic approximations, depending on the product and formulation.

The exact formula matters less than the verification principle:

increase v_max and check whether the target price changes materially.

If it does, the original variance domain was too small or the boundary rule was inappropriate.

14. The far asset-price boundary must also be tested

For a vanilla call, large-S asymptotics can guide the boundary. Barrier and digital payoffs can require different treatment.

As in the one-dimensional Black–Scholes PDE, domain truncation is a source of model-independent numerical error. A sophisticated stochastic-volatility model does not excuse a poorly chosen computational boundary.

15. Payoff kinks can create temporal oscillations

A vanilla call payoff has a kink at the strike. Barrier and digital contracts can be even less smooth.

Second-order schemes can react badly to such nonsmooth terminal data at the first few time steps. Oscillations may appear in price derivatives even when the option value itself looks plausible.

Rannacher smoothing is a common remedy: begin with a few strongly damping backward-Euler half-steps, then switch to the higher-order ADI scheme.

The purpose is not cosmetic. Greeks such as delta and gamma can be more sensitive to high-frequency numerical error than the price itself.

16. Mixed-derivative discretisation is a weak link

On a rectangular grid, a central approximation to V_Sv uses values from diagonal neighbours. On non-uniform grids, the coefficients become more complicated.

Errors in the mixed derivative are easy to hide because the solver can still converge numerically—to the wrong discrete equation.

A useful unit test sets ρ = 0. The mixed derivative should vanish exactly from the assembled operator. Then restore nonzero ρ and compare against a trusted Heston transform price for European options.

17. Inputs and outputs

Inputs can include:

• spot S₀;
• initial variance v₀;
• κ, θ, σ_v and ρ;
• discount and dividend/foreign-rate curves;
• payoff and exercise style;
• S and v grid transformations;
• domain maxima;
• number of spatial nodes;
• time-step count;
• ADI scheme and θ;
• boundary conditions;
• Rannacher startup policy;
• interpolation rule.

Outputs can include:

• price at the requested state;
• full S–v value surface;
• delta, gamma and variance sensitivity;
• exercise or barrier diagnostics;
• temporal and spatial convergence tables;
• boundary-sensitivity results;
• ADI scheme comparison;
• linear-solve residuals;
• runtime and memory use.

18. Benchmark against the Heston transform formula

For European vanilla options, Heston’s characteristic-function solution provides a powerful reference.

The PDE engine and transform engine use the same model parameters but very different numerical machinery. Agreement across strikes, maturities, correlations and Feller regimes is strong evidence that both implementations are behaving sensibly.

One matching price is not enough. A bug can cancel accidentally at one point.

19. Convergence should be measured in two spatial directions and time

A two-dimensional grid has at least three independent resolution controls:

• ΔS or its non-uniform analogue;
• Δv;
• Δτ.

Refining all three simultaneously can hide which source dominates the error. A better diagnostic changes one dimension at a time and then performs combined refinement.

Convergence should also be tested for Greeks, not only prices.

20. Evidence polarity

Evidence for confidence includes close agreement with transform prices for European vanillas, smooth convergence under independent refinement of price, variance and time grids, small sensitivity to farther domain boundaries, stable results in both Feller-satisfied and Feller-violated cases, consistent output across MCS/HV at fine grids, and sensible no-arbitrage shape behaviour.

Evidence against confidence includes prices that depend strongly on the artificial variance ceiling, oscillatory Greeks after maturity, failure when |ρ| becomes large, scheme-dependent answers that do not converge together, discontinuous behaviour as v₀ approaches zero, negative vanilla option values, mixed-derivative errors revealed by ρ = 0 tests, or apparent stability without grid convergence.

21. Counterexample: stable but inaccurate

An ADI scheme can remain bounded on a coarse grid and therefore appear “stable.” Yet the price may still be materially wrong because the strike region has too few nodes or the variance boundary is too close.

Falsifier: double the spatial resolution while holding the time scheme fixed. If the result moves materially, stability was not accuracy.

22. Counterexample: correlation exposes a hidden implementation bug

A solver may match benchmark prices when ρ = 0 but fail badly at ρ = −0.8.

The likely weak links include the mixed-derivative stencil, operator split or correlation sign convention.

Falsifier: sweep ρ from negative to positive values and compare every point with an independent transform pricer.

23. Counterexample: Feller violation changes the error regime

A solver calibrated and tested only where 2κθ is much larger than σ_v² may look excellent. A market calibration that pushes the variance process toward zero can expose boundary and grid weaknesses.

Falsifier: run parameter sets on both sides of the Feller boundary and refine the low-variance grid independently.

24. Counterexample: a perfect PDE solution can still be a poor market model

Suppose the ADI price agrees with the transform solution to six decimal places. That proves the numerical PDE solver is solving the Heston model accurately.

It does not prove Heston fits the observed volatility surface, forward skew dynamics or hedging behaviour.

Falsifier: separate numerical error from calibration error. If the PDE and transform agree but market residuals remain structured, the model—not the solver—is the weak link.

25. Alternatives

Characteristic-function integration is often faster for European vanillas under Heston.

Monte Carlo is flexible for path dependence and higher dimensions, but requires simulation error control and careful variance-process discretisation.

Fourier FFT/COS methods exploit characteristic functions and can price grids of strikes efficiently.

Finite elements provide flexible meshes and alternative convergence properties.

ADI PDE methods are especially attractive when the state dimension is low and contract features are naturally represented as boundaries or exercise constraints.

26. Connections to the surrounding Bukit Timah Tutor knowledge estate

The one-dimensional numerical foundation is finite-difference Black–Scholes PDE algorithms. This page extends that logic to a genuinely two-dimensional stochastic-volatility problem with a mixed derivative.

Parameter estimation belongs to Heston calibration algorithms. The PDE solver should consume calibrated risk-neutral parameters rather than silently refit them.

Simulation provides an independent numerical route. Andersen QE algorithms address variance simulation, while Monte Carlo pricing provides the broader path engine.

Transform alternatives include Carr–Madan FFT and Fourier-COS pricing.

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

27. What would falsify confidence?

Confidence should be withdrawn if benchmark vanilla prices do not converge to an independent Heston transform implementation; if the mixed term fails the ρ = 0 test; if results depend materially on arbitrary domain ceilings; if Greeks oscillate under reasonable refinement; if Feller-violated calibrations destabilise the grid; if different second-order ADI schemes converge to different values; or if numerical agreement is mistaken for empirical validity.

28. Verification and update triggers

Preserve the Heston parameters, curves, payoff, grid maps, domain limits, ADI scheme, θ value, startup smoothing, boundary rules and every convergence table with each model version.

Revalidate after changes to grid construction, mixed-derivative stencils, boundary logic, linear solvers, Heston parameter ranges, exercise features, barrier handling or numerical libraries. Trigger review when calibration begins producing more extreme correlation, volatility-of-variance or near-zero-variance regimes than the original test suite covered.

29. Primary and high-quality references

• Steven L. Heston, A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options, Review of Financial Studies, 1993.
• Karel J. in ’t Hout and Sven Foulon, ADI Finite Difference Schemes for Option Pricing in the Heston Model with Correlation, International Journal of Numerical Analysis and Modeling, 2010.
• Fabrice Douglas Rouah, Methods for Finite Differences, in The Heston Model and Its Extensions in Matlab and C#, Wiley, 2013.
• Tinne Haentjens and Karel J. in ’t Hout, ADI Schemes for Pricing American Options under the Heston Model, Applied Mathematical Finance, 2015.
• Bénédicte Alziary and Peter Takáč, The Heston stochastic volatility model has a boundary trace at zero volatility, 2023.

Educational boundary: ADI is numerical machinery for solving a chosen stochastic-volatility PDE. Accurate numerics are necessary for trustworthy pricing, but they do not make the model assumptions true.