Robertson stiff chemical kinetics simulation

What this research found

The Robertson reaction network — three chemical species whose rate constants span nine orders of magnitude — is the standard benchmark for stiff differential equations. K-Dense integrated it from time zero to 100,000 with an adaptive implicit backward differentiation formula solver, then built explicit integrators from scratch to measure exactly how badly they fail. The implicit method crossed the whole interval in 377 accepted steps while conserving mass to 1.33 x 10^-15; an explicit method would need roughly 350 to 490 million steps to do the same.

• The implicit solver reached the end of the interval in 377 accepted steps, using 1,068 right-hand-side evaluations, 13 Jacobian evaluations, and 81 LU factorizations. Its step size grew by more than eleven orders of magnitude, from about 6.6 x 10^-8 during the initial transient to about 3.6 x 10^3 on the slow manifold.
• The slow reactant falls from 1 to 0.96646 at time 1, 0.61724 at 100, 0.10730 at 10,000, and 0.017866 at 100,000, while the terminal product rises correspondingly. The reactive intermediate never exceeds about 3.6 x 10^-5 and decays to roughly 7 x 10^-8.
• Total concentration, an invariant never imposed on the solver, held to a maximum absolute deviation of 1.33 x 10^-15 with a mean of 7.33 x 10^-16 — the level of double-precision unit round-off, indicating the error is dominated by floating-point round-off rather than integration error.
• Contrary to the common assumption that stiff problems are hardest at the start, the Jacobian spectrum at time zero is {-0.04, 0, 0} and the system is not stiff there at all. The dominant eigenvalue modulus climbs to about 2.19 x 10^3 at the intermediate's peak and to a global maximum of about 9.83 x 10^3 near the end of the run, giving an effective stiffness ratio exceeding 2 x 10^5.
• Measured stability limits confirm the eigenvalue theory. Forward Euler stayed stable up to a step of 7.70 x 10^-4 and classical fourth-order Runge-Kutta up to 1.16 x 10^-3, matching the predicted bounds to within 12% and 21%. Just past those thresholds both blow up abruptly, forward Euler diverging near time 25.2 and Runge-Kutta near 31.6.
• Adaptive explicit solvers do not escape the problem. Dormand-Prince and DOP853 reached time 1 easily in a few thousand evaluations but both hit a 30-second wall-clock guard before reaching 100,000, because adaptivity controls accuracy and cannot lift the stability ceiling.

How it was done

The three-species system was integrated with a variable-order, variable-step backward differentiation formula method supplied with the Jacobian in closed form, using a relative tolerance of 10^-6 and per-component absolute tolerances of 10^-6, 10^-14, and 10^-6. The tight middle value places the error floor about eight orders of magnitude below the trace intermediate's peak, so that species is controlled by relative rather than absolute error and the fast quasi-steady balance stays accurate. Stiffness was then quantified three ways: by evaluating the Jacobian spectrum along the reference trajectory, by bisection search for the largest stable fixed step of hand-written forward-Euler and Runge-Kutta integrators, and by direct blow-up demonstrations just above and below that threshold. Cost to reach the end of the interval with an explicit method was extrapolated from measured per-step timings, and the log-spaced trajectory plus the full adaptive step history were exported so every figure can be regenerated.

Data sources

• Robertson, in Numerical Analysis: An Introduction, Academic Press (1966) — the original reaction-rate system
• Hairer & Wanner, Solving Ordinary Differential Equations II, 2nd edition (1996) — stiff and differential-algebraic problems
• Gear, Numerical Initial Value Problems in Ordinary Differential Equations (1971) — backward differentiation formula methods
• Enright, Hull & Lindberg, BIT Numerical Mathematics 15:10 (1975) — stiff-solver comparisons