Notebooks · lib/ode/stiff_reference

lib/ode against the stiff test set's published reference solutions

Five stiff problems, each solved by Radau IIA and by BDF and compared with the reference solution published with the problem: the IVP test set of Mazzia and Magherini (Bari, release 2.4, 2008) for Robertson to t = 10¹¹, Van der Pol with eps = 10⁻⁶ (mu = 1000), HIRES and the Oregonator (whose values agree with those in Hairer and Wanner, Solving ODEs II, IV.10). The error is the largest relative error over the components at the end time; the bound is 100 rtol for Radau and 1000 rtol for BDF, whose global error on the Oregonator grows over its oscillations (SciPy's BDF makes the same 1.5 × 10⁻⁶ at rtol = 10⁻⁸).

Van der Pol with eps = 10⁻¹² (mu = 10⁶) has no published value; Radau's result is held to SciPy's Radau at rtol = 10⁻¹⁰ (1.705546217535324, python/tests/test_ode.py runs both). BDF fails on it at t = 0.807, at the first relaxation jump, as SciPy's BDF does: the step size it needs there underflows.

E5 is not here: at its end time y2 and y3 are 10⁻²² and their difference 10⁻²⁵, and no setting of rtol and atol makes SciPy's or these Radau and BDF reproduce the reference (see the report on branch ode).

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 26–105

from ode.solve import ode_cfg, ode_solve, no_events
from ode.core import BDF, DONE, RADAU, JAC_GIVEN
from ode.problems import (
    robertson, robertson_jac, vdpol, vdpol_jac, hires, hires_jac, oregonator, oregonator_jac,
)

RTOL = 1.0e-8
# the error a test allows, in units of rtol: BDF's global error on the
# Oregonator grows over its oscillations (SciPy's BDF: 1.5e-6 at rtol 1e-8)
ALLOWED = 100.0
ALLOWED_BDF = 1000.0
MAX_STEPS: Int = 1000000
NO_PARAMS = zeros(0)

ROBER_T = 1.0e11
ROBER_Y0 = [1.0, 0.0, 0.0]
ROBER_P = [0.04, 1.0e4, 3.0e7]
ROBER_REF = [0.2083340149701255e-07, 0.8333360770334713e-13, 0.9999999791665050]
VDPOL_T = 2.0
VDPOL_Y0 = [2.0, 0.0]
VDPOL_REF = [0.1706167732170483e+01, -0.8928097010247975]
HIRES_T = 321.8122
HIRES_Y0 = [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0057]
HIRES_REF = [0.7371312573325668e-3, 0.1442485726316185e-3, 0.5888729740967575e-4, 0.1175651343283149e-2,
             0.2386356198831331e-2, 0.6238968252742796e-2, 0.2849998395185769e-2, 0.2850001604814231e-2]
OREGO_T = 360.0
OREGO_Y0 = [1.0, 2.0, 3.0]
OREGO_REF = [0.1000814870318523e+01, 0.1228178521549917e+04, 0.1320554942846706e+03]
VDP6_REF = [1.705546217535324, -0.8934763625465864]

#| The largest relative error of y against ref.
def rel_err(y: Arr, ref: Arr) -> Num {
    return amax(abs(y - ref) / abs(ref))
}

#| Solve to t1 and return the state there, the status and the step count.
def end_state(f, jac, method: Int, t1: Num, y0: Arr, p: Arr, rtol: Num, atol: Num) {
    c = ode_cfg(method, t1, rtol, atol * ones(len(y0)), 0.0, inf, MAX_STEPS, JAC_GIVEN)
    r = ode_solve(f, jac, no_events, 0.0, t1, y0, p, c, [t1], zeros(0), zeros(0))
    return r.y[0, :], r.status, r.nsteps
}

methods = [RADAU, BDF]
# Robertson, Van der Pol, HIRES, the Oregonator
NPROBLEMS: Int = 4
errs = zeros(len(methods), NPROBLEMS)
done = 0
mi = 0
for mm in methods {
    m: Int = int(mm)
    y1, st1, n1 = end_state(robertson, robertson_jac, m, ROBER_T, ROBER_Y0, ROBER_P, RTOL, 1.0e-20)
    y2, st2, n2 = end_state(vdpol, vdpol_jac, m, VDPOL_T, VDPOL_Y0, [1.0e-6], RTOL, 1.0e-8)
    y3, st3, n3 = end_state(hires, hires_jac, m, HIRES_T, HIRES_Y0, NO_PARAMS, RTOL, 1.0e-12)
    y4, st4, n4 = end_state(oregonator, oregonator_jac, m, OREGO_T, OREGO_Y0, NO_PARAMS, RTOL, 1.0e-8)
    errs[mi, :] = [rel_err(y1, ROBER_REF), rel_err(y2, VDPOL_REF), rel_err(y3, HIRES_REF), rel_err(y4, OREGO_REF)]
    print(f"method {m}: steps {n1} {n2} {n3} {n4}; relative errors {errs[mi, 0]:.2e} {errs[mi, 1]:.2e} {errs[mi, 2]:.2e} {errs[mi, 3]:.2e}")
    if st1 == DONE and st2 == DONE and st3 == DONE and st4 == DONE {
        done += 1
    }
    mi += 1
}
y6, st6, n6 = end_state(vdpol, vdpol_jac, RADAU, VDPOL_T, VDPOL_Y0, [1.0e-12], RTOL, 1.0e-8)
e6 = rel_err(y6, VDP6_REF)
print(f"Van der Pol mu = 1e6, Radau: {n6} steps, relative difference from SciPy {e6:.2e}")

test "Radau IIA and BDF reach the end of all four problems" {
    assert_near(done, 2.0, 0.0)
}
test "Radau IIA matches the published reference solutions to 100 rtol" {
    assert_lt(amax(errs[0, :]), ALLOWED * RTOL)
}
test "BDF matches the published reference solutions to 1000 rtol" {
    assert_lt(amax(errs[1, :]), ALLOWED_BDF * RTOL)
}
test "Radau IIA solves Van der Pol at mu = 1e6 as SciPy's Radau does" {
    assert_lt(e6, ALLOWED * RTOL)
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

from ode.solve import ode_cfg, ode_solve, no_events
from ode.core import BDF, DONE, RADAU, JAC_GIVEN
from ode.problems import (
    robertson, robertson_jac, vdpol, vdpol_jac, hires, hires_jac, oregonator, oregonator_jac,
)

RTOL = 1.0e-8
# the error a test allows, in units of rtol: BDF's global error on the
# Oregonator grows over its oscillations (SciPy's BDF: 1.5e-6 at rtol 1e-8)
ALLOWED = 100.0
ALLOWED_BDF = 1000.0
MAX_STEPS: Int = 1000000
NO_PARAMS = zeros(0)

ROBER_T = 1.0e11
ROBER_Y0 = [1.0, 0.0, 0.0]
ROBER_P = [0.04, 1.0e4, 3.0e7]
ROBER_REF = [0.2083340149701255e-07, 0.8333360770334713e-13, 0.9999999791665050]
VDPOL_T = 2.0
VDPOL_Y0 = [2.0, 0.0]
VDPOL_REF = [0.1706167732170483e+01, -0.8928097010247975]
HIRES_T = 321.8122
HIRES_Y0 = [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0057]
HIRES_REF = [0.7371312573325668e-3, 0.1442485726316185e-3, 0.5888729740967575e-4, 0.1175651343283149e-2,
             0.2386356198831331e-2, 0.6238968252742796e-2, 0.2849998395185769e-2, 0.2850001604814231e-2]
OREGO_T = 360.0
OREGO_Y0 = [1.0, 2.0, 3.0]
OREGO_REF = [0.1000814870318523e+01, 0.1228178521549917e+04, 0.1320554942846706e+03]
VDP6_REF = [1.705546217535324, -0.8934763625465864]

#| The largest relative error of y against ref.
def rel_err(y: Arr, ref: Arr) -> Num {
    return amax(abs(y - ref) / abs(ref))
}

#| Solve to t1 and return the state there, the status and the step count.
def end_state(f, jac, method: Int, t1: Num, y0: Arr, p: Arr, rtol: Num, atol: Num) {
    c = ode_cfg(method, t1, rtol, atol * ones(len(y0)), 0.0, inf, MAX_STEPS, JAC_GIVEN)
    r = ode_solve(f, jac, no_events, 0.0, t1, y0, p, c, [t1], zeros(0), zeros(0))
    return r.y[0, :], r.status, r.nsteps
}

methods = [RADAU, BDF]
# Robertson, Van der Pol, HIRES, the Oregonator
NPROBLEMS: Int = 4
errs = zeros(len(methods), NPROBLEMS)
done = 0
mi = 0
for mm in methods {
    m: Int = int(mm)
    y1, st1, n1 = end_state(robertson, robertson_jac, m, ROBER_T, ROBER_Y0, ROBER_P, RTOL, 1.0e-20)
    y2, st2, n2 = end_state(vdpol, vdpol_jac, m, VDPOL_T, VDPOL_Y0, [1.0e-6], RTOL, 1.0e-8)
    y3, st3, n3 = end_state(hires, hires_jac, m, HIRES_T, HIRES_Y0, NO_PARAMS, RTOL, 1.0e-12)
    y4, st4, n4 = end_state(oregonator, oregonator_jac, m, OREGO_T, OREGO_Y0, NO_PARAMS, RTOL, 1.0e-8)
    errs[mi, :] = [rel_err(y1, ROBER_REF), rel_err(y2, VDPOL_REF), rel_err(y3, HIRES_REF), rel_err(y4, OREGO_REF)]
    print(f"method {m}: steps {n1} {n2} {n3} {n4}; relative errors {errs[mi, 0]:.2e} {errs[mi, 1]:.2e} {errs[mi, 2]:.2e} {errs[mi, 3]:.2e}")
    if st1 == DONE and st2 == DONE and st3 == DONE and st4 == DONE {
        done += 1
    }
    mi += 1
}
y6, st6, n6 = end_state(vdpol, vdpol_jac, RADAU, VDPOL_T, VDPOL_Y0, [1.0e-12], RTOL, 1.0e-8)
e6 = rel_err(y6, VDP6_REF)
print(f"Van der Pol mu = 1e6, Radau: {n6} steps, relative difference from SciPy {e6:.2e}")

test "Radau IIA and BDF reach the end of all four problems" {
    assert_near(done, 2.0, 0.0)
}
test "Radau IIA matches the published reference solutions to 100 rtol" {
    assert_lt(amax(errs[0, :]), ALLOWED * RTOL)
}
test "BDF matches the published reference solutions to 1000 rtol" {
    assert_lt(amax(errs[1, :]), ALLOWED_BDF * RTOL)
}
test "Radau IIA solves Van der Pol at mu = 1e6 as SciPy's Radau does" {
    assert_lt(e6, ALLOWED * RTOL)
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors