Notebooks · lib/ode/gallery_robertson

Robertson's stiff kinetics, drawn

The first figure of vipyr's gallery example plot_ode_robertson.py, drawn in Viv: Robertson's three-species reaction (rates 0.04, 10⁴ and 3 × 10⁷) solved by Radau IIA from t = 0 to 10¹¹ at rtol 10⁻⁸, through ode_solve in lib/ode/solve.viv with the Jacobian from grad. The middle species is drawn times 1e4. The end state is checked against the reference solution the IVP test set of Mazzia and Magherini (Bari) publishes, as vipyr's script checks it.

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 13–58

from ode.core import RADAU, JAC_GIVEN
from ode.solve import ode_cfg, ode_solve, no_events
from ode.problems import robertson, robertson_jac

RATES = [0.04, 1e4, 3e7]
T_END: Num = 1e11
RTOL: Num = 1e-8
ATOL: Num = 1e-20
#: output at 161 times, log-spaced over 16 decades
N_OUT: Int = 161
LOG_T_FIRST: Num = -5.0
LOG_T_LAST: Num = 11.0
#: the middle species is drawn scaled up, as vipyr draws it
Y2_SCALE: Num = 1e4
#: the Bari IVP test set's reference at t = 1e11
Y_REF = [0.2083340149701255e-07, 0.8333360770334713e-13, 0.9999999791665050]
REF_REL: Num = 1e-4
#: the first step chosen by the solver, and no cap on the step
AUTO_FIRST: Num = 0.0
MAX_STEPS: Int = 100000

t_out = [10.0 ^ x for x in linspace(LOG_T_FIRST, LOG_T_LAST, N_OUT)]
c = ode_cfg(RADAU, T_END, RTOL, ones(3) * ATOL, AUTO_FIRST, inf, MAX_STEPS, JAC_GIVEN)
r = ode_solve(robertson, robertson_jac, no_events, 0.0, T_END, [1.0, 0.0, 0.0], RATES, c, t_out, zeros(0),
              zeros(0))
show r.nsteps
last = r.y[len(t_out) - 1, :]
show last

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

test "the end state is the reference solution's" {
    for i to len(Y_REF) {
        assert_lt(abs(last[i] - Y_REF[i]) / Y_REF[i], REF_REL)
    }
}

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

from ode.core import RADAU, JAC_GIVEN
from ode.solve import ode_cfg, ode_solve, no_events
from ode.problems import robertson, robertson_jac

RATES = [0.04, 1e4, 3e7]
T_END: Num = 1e11
RTOL: Num = 1e-8
ATOL: Num = 1e-20
#: output at 161 times, log-spaced over 16 decades
N_OUT: Int = 161
LOG_T_FIRST: Num = -5.0
LOG_T_LAST: Num = 11.0
#: the middle species is drawn scaled up, as vipyr draws it
Y2_SCALE: Num = 1e4
#: the Bari IVP test set's reference at t = 1e11
Y_REF = [0.2083340149701255e-07, 0.8333360770334713e-13, 0.9999999791665050]
REF_REL: Num = 1e-4
#: the first step chosen by the solver, and no cap on the step
AUTO_FIRST: Num = 0.0
MAX_STEPS: Int = 100000

t_out = [10.0 ^ x for x in linspace(LOG_T_FIRST, LOG_T_LAST, N_OUT)]
c = ode_cfg(RADAU, T_END, RTOL, ones(3) * ATOL, AUTO_FIRST, inf, MAX_STEPS, JAC_GIVEN)
r = ode_solve(robertson, robertson_jac, no_events, 0.0, T_END, [1.0, 0.0, 0.0], RATES, c, t_out, zeros(0),
              zeros(0))
show r.nsteps
last = r.y[len(t_out) - 1, :]
show last

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

test "the end state is the reference solution's" {
    for i to len(Y_REF) {
        assert_lt(abs(last[i] - Y_REF[i]) / Y_REF[i], REF_REL)
    }
}

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