Stiff Systems with CVODE BDF¶
Stiff ODEs contain dynamics on vastly different time scales. Explicit methods (Adams) need tiny time steps to remain stable, while implicit methods (BDF) can take much larger steps.
CVODE's BDF method with dense direct linear solver is well-suited for moderate-dimensional stiff systems.
load("numerics")$
load("numerics-sundials")$
load("ax-plots")$
1. Van der Pol Oscillator¶
The Van der Pol oscillator becomes increasingly stiff as $\mu$ grows: $$\ddot{x} - \mu(1 - x^2)\dot{x} + x = 0$$
Written as a system: $\dot{x} = y$, $\dot{y} = \mu(1-x^2)y - x$.
/* Moderate stiffness: mu = 10 */
mu : 10$
tspan1 : np_linspace(0.0, 60.0, 5000)$
res1 : np_cvode([y, mu*(1 - x^2)*y - x],
[t, x, y], [2.0, 0.0], tspan1,
1e-8, 1e-8, bdf)$
t1 : np_slice(res1, all, 0)$
x1 : np_slice(res1, all, 1)$
y1 : np_slice(res1, all, 2)$
ax_draw2d(
color="#2563eb", line_width=1.5,
lines(t1, x1),
title=printf(false, "Van der Pol Oscillator (mu=~d)", mu),
xlabel="t", ylabel="x",
grid=true
)$
/* Phase portrait shows the relaxation oscillation limit cycle */
ax_draw2d(
color="#7c3aed", line_width=1,
lines(x1, y1),
title="Van der Pol Phase Portrait",
xlabel="x", ylabel="dx/dt",
grid=true
)$
Highly stiff: $\mu = 1000$¶
With $\mu = 1000$ the system has sharp transitions where $x$ jumps rapidly between $\pm 2$. BDF handles this where Adams would struggle.
mu2 : 1000$
tspan2 : np_linspace(0.0, 6000.0, 20000)$
res2 : np_cvode([y2, mu2*(1 - x2^2)*y2 - x2],
[t, x2, y2], [2.0, 0.0], tspan2,
1e-6, 1e-9, bdf)$
t2 : np_slice(res2, all, 0)$
x2_col : np_slice(res2, all, 1)$
ax_draw2d(
color="#dc2626", line_width=1,
lines(t2, x2_col),
title=printf(false, "Van der Pol (mu=~d) — Relaxation Oscillation", mu2),
xlabel="t", ylabel="x",
grid=true
)$
2. Robertson Chemical Kinetics¶
A classic stiff test problem (SUNDIALS cvRoberts_dns):
$$y_1' = -0.04\,y_1 + 10^4\,y_2\,y_3$$
$$y_2' = 0.04\,y_1 - 10^4\,y_2\,y_3 - 3{\times}10^7\,y_2^2$$
$$y_3' = 3{\times}10^7\,y_2^2$$
The three species evolve on time scales spanning many orders of magnitude.
/* Logarithmically spaced output times from 0.001 to 1e6 */
tspan3 : append([0.0],
makelist(10^(k/10.0), k, -30, 60))$
res3 : np_cvode(
[-0.04*s1 + 1e4*s2*s3,
0.04*s1 - 1e4*s2*s3 - 3e7*s2^2,
3e7*s2^2],
[t, s1, s2, s3], [1.0, 0.0, 0.0], tspan3,
1e-8, 1e-12, bdf)$
t3 : np_slice(res3, all, 0)$
s1_col : np_slice(res3, all, 1)$
s2_col : np_slice(res3, all, 2)$
s3_col : np_slice(res3, all, 3)$
n3 : first(np_shape(res3)) - 1$
printf(true, "Final state at t=~,0e:~%", last(tspan3))$
printf(true, " y1 = ~,6e~%", np_ref(res3, n3, 1))$
printf(true, " y2 = ~,6e~%", np_ref(res3, n3, 2))$
printf(true, " y3 = ~,6e~%", np_ref(res3, n3, 3))$
Final state at t=1.e+6:
o143 y1 = 2.031484e-3
o143 y2 = 8.142279e-9
o143 y3 = 9.979685e-1
o143$$\mathbf{false}$$
false
/* Plot on log10(t) scale — skip t=0 to avoid log(0) */
nr3 : first(np_shape(res3))$
log_t3 : np_log10(np_slice(res3, [1, nr3], 0))$
s1_skip : np_slice(res3, [1, nr3], 1)$
s3_skip : np_slice(res3, [1, nr3], 3)$
ax_draw2d(
color="#2563eb", line_width=2, name="y1",
lines(log_t3, s1_skip),
color="#16a34a", line_width=2, name="y3",
lines(log_t3, s3_skip),
title="Robertson Kinetics",
xlabel="log10(t)", ylabel="concentration",
grid=true, showlegend=true
)$
/* y2 is tiny — scale by 1e4 for visibility */
s2_skip : np_slice(res3, [1, nr3], 2)$
y2_scaled : np_scale(1e4, s2_skip)$
ax_draw2d(
color="#f59e0b", line_width=2, name="y2 x 10^4",
lines(log_t3, y2_scaled),
title="Robertson Kinetics — y2 (scaled)",
xlabel="log10(t)", ylabel="y2 x 10^4",
grid=true, showlegend=true
)$
3. ROBER Consistency Check¶
The Robertson system conserves total mass: $y_1 + y_2 + y_3 = 1$ for all $t$.
mass_sum : np_add(s1_col, np_add(s2_col, s3_col))$
nr : first(np_shape(res3))$
ones : np_ones([nr])$
printf(true, "Max deviation from mass=1: ~,2e~%",
np_max(np_abs(np_sub(mass_sum, ones))))$
Max deviation from mass=1: 1.66e-13
o171$$\mathbf{false}$$
false