Robertson's stiff kinetics over sixteen decades

Three species, rate constants 0.04, 10⁴ and 3 × 10⁷: the fast reactions settle in microseconds and the slow one keeps going until t = 1e11. An explicit method would need steps of about 10⁻⁴ the whole way. vp.ode.solve with Radau IIA and with BDF crosses the range in under two thousand steps, lands on the published reference solution, and takes the steps SciPy's solve_ivp takes.

The problem

"robertson" is a def of lib/ode/problems.viv:

y1' = -k1 y1 + k2 y2 y3
y2' =  k1 y1 - k2 y2 y3 - k3 y2^2
y3' =  k3 y2^2

The Jacobian is exact: the driver writes one reverse-mode grad per row. Output at 162 times, t = 0 and ten per decade from 10⁻⁵ to 10¹¹, by the dense output of each step.

import numpy as np
import vipyr as vp
from vipyr.figure import figure

k = [0.04, 1e4, 3e7]
t_out = np.concatenate([[0.0], np.logspace(-5, 11, 161)])
runs = {}
for method in ("radau", "bdf"):
    runs[method] = vp.ode.solve(
        "robertson",
        (0.0, 1e11),
        [1.0, 0.0, 0.0],
        method=method,
        rtol=1e-8,
        atol=1e-20,
        params=k,
        t_eval=t_out,
    )
    s = runs[method].stats
    print(
        "%-6s %5d steps, %5d rhs evaluations, %4d Jacobians, %4d LU "
        "factorizations, %3d rejected"
        % (method, s["steps"], s["nfev"], s["njev"], s["nlu"], s["rejected"])
    )

radau   1717 steps, 13107 rhs evaluations,  357 Jacobians,  988 LU factorizations,   2 rejected
bdf     1757 steps,  5411 rhs evaluations,   25 Jacobians,  324 LU factorizations,   2 rejected

Against the reference solution

The IVP test set of Mazzia and Magherini (Bari) publishes the state at t = 10¹¹ to sixteen digits. At rtol = 10⁻⁸ Radau IIA is within 10⁻¹⁰ of it relative to each component, and BDF within 2e-7. Mass is conserved: y1 + y2 + y3 = 1 throughout.

ref = np.array(
    [0.2083340149701255e-07, 0.8333360770334713e-13, 0.9999999791665050]
)
for method, r in runs.items():
    err = np.abs(r.y[:, -1] - ref) / ref
    print(
        "%-6s y(1e11) = %s   largest relative error %.1e   max |sum - 1| = "
        "%.1e"
        % (
            method,
            np.array2string(r.y[:, -1], precision=10),
            err.max(),
            np.abs(r.y.sum(axis=0) - 1).max(),
        )
    )

radau  y(1e11) = [2.0833401496e-08 8.3333607699e-14 9.9999997917e-01]   largest relative error 5.4e-11   max |sum - 1| = 1.2e-15
bdf    y(1e11) = [2.0833404173e-08 8.3333618409e-14 9.9999997917e-01]   largest relative error 1.3e-07   max |sum - 1| = 2.0e-15

The same steps as SciPy

lib/ode/radau.viv and lib/ode/bdf.viv follow SciPy's step-size, Newton and order logic, so they take SciPy's steps (to within a handful) and agree with its answer at the end to rounding.

from scipy.integrate import solve_ivp


def rob(t, y):
    return [
        -k[0] * y[0] + k[1] * y[1] * y[2],
        k[0] * y[0] - k[1] * y[1] * y[2] - k[2] * y[1] ** 2,
        k[2] * y[1] ** 2,
    ]


def rob_jac(t, y):
    return [
        [-k[0], k[1] * y[2], k[1] * y[1]],
        [k[0], -k[1] * y[2] - 2 * k[2] * y[1], -k[1] * y[1]],
        [0.0, 2 * k[2] * y[1], 0.0],
    ]


for method, name in (("radau", "Radau"), ("bdf", "BDF")):
    sp = solve_ivp(
        rob,
        (0.0, 1e11),
        [1.0, 0.0, 0.0],
        method=name,
        rtol=1e-8,
        atol=1e-20,
        jac=rob_jac,
    )
    d = np.abs(runs[method].y[:, -1] - sp.y[:, -1]) / sp.y[:, -1]
    print(
        "%-6s SciPy %5d steps, vp.ode %5d steps; end states differ by %.0e "
        "relative"
        % (method, len(sp.t) - 1, runs[method].stats["steps"], d.max())
    )

radau  SciPy  1718 steps, vp.ode  1717 steps; end states differ by 9e-13 relative
bdf    SciPy  1757 steps, vp.ode  1757 steps; end states differ by 1e-13 relative

The picture

The species against log time, y2 scaled by 10⁴ as is customary. The intermediate y2 jumps to its quasi-steady level within 10⁻⁴ time units, then decays with y1 for ten decades while y3 takes over.

r = runs["radau"]
mask = t_out > 0
fig = figure(
    "Robertson's kinetics (Radau IIA, rtol 1e-8)", "t", "concentration"
)
fig.xscale("log")
fig.plot(t_out[mask], r.y[0, mask], "y1")
fig.plot(t_out[mask], 1e4 * r.y[1, mask], "1e4 y2")
fig.plot(t_out[mask], r.y[2, mask], "y3")
fig.style("y1", "s1 nomarker")
fig.style("1e4 y2", "s2 nomarker")
fig.style("y3", "s3 nomarker")

Source

gallery/math/plot_ode_robertson.py, run as written when these pages were built. The output and the figures above are from that run.

The Viv program it calls is lib/ode/gallery_robertson.viv: its atlas, the defs and the calls between them, or open it in the workbench.