ODE Integration with SUNDIALS CVODE¶
np_cvode wraps the SUNDIALS CVODE solver for initial value problems $\dot{y} = f(t, y)$.
Two calling conventions:
- Expression mode:
np_cvode([f1, ...], [t, y1, ...], y0, tspan) - Function mode:
np_cvode(f, y0, tspan)wheref(t, y)returns a list
Supports Adams (non-stiff, default) and BDF (stiff) methods.
load("numerics")$
load("numerics-sundials")$
load("ax-plots")$
1. Exponential Decay¶
$\dot{y} = -y$, $y(0) = 1$. Exact solution: $y(t) = e^{-t}$.
tspan : np_linspace(0.0, 5.0, 100)$
result : np_cvode([-y], [t, y], [1.0], tspan)$
t_col : np_slice(result, all, 0)$
y_col : np_slice(result, all, 1)$
y_exact : np_exp(np_scale(-1.0, t_col))$
print("Max error:", np_max(np_abs(np_sub(y_col, y_exact))))$
Max error: 1.3237647200625702e-8
ax_draw2d(
color="#2563eb", line_width=2, name="CVODE",
lines(t_col, y_col),
color="#dc2626", dash="dash", name="exp(-t)",
explicit(exp(-t), t, 0, 5),
title="Exponential Decay", xlabel="t", ylabel="y",
grid=true, showlegend=true
)$
2. Harmonic Oscillator¶
$\dot{x} = v$, $\dot{v} = -x$, with $(x_0, v_0) = (1, 0)$.
Exact: $x(t) = \cos t$, $v(t) = -\sin t$. Phase portrait should be a perfect circle.
tspan2 : np_linspace(0.0, float(4*%pi), 400)$
res2 : np_cvode([v, -x], [t, x, v], [1.0, 0.0], tspan2)$
t2 : np_slice(res2, all, 0)$
x2 : np_slice(res2, all, 1)$
v2 : np_slice(res2, all, 2)$
ax_draw2d(
color="#2563eb", line_width=2, name="x(t)",
lines(t2, x2),
color="#f59e0b", line_width=2, name="v(t)",
lines(t2, v2),
title="Harmonic Oscillator", xlabel="t", ylabel="",
grid=true, showlegend=true
)$
/* Phase portrait */
ax_draw2d(
color="#7c3aed", line_width=2,
lines(x2, v2),
title="Phase Portrait", xlabel="x", ylabel="v",
grid=true
)$
3. Function Mode¶
The same oscillator using a callback function instead of symbolic expressions.
osc_rhs(t, y) := [np_ref(y, 1), -np_ref(y, 0)]$
res3 : np_cvode(osc_rhs, [1.0, 0.0], np_linspace(0.0, float(2*%pi), 200))$
np_shape(res3);
4. Stiff System — Adams vs BDF¶
$\dot{y} = -1000\,y + 3000 - 2000\,e^{-t}$, $y(0) = 0$.
Exact: $y(t) = 3 - 0.998\,e^{-1000t} - 2.002\,e^{-t}$.
The fast transient ($e^{-1000t}$) decays in ~0.01 time units, making this stiff. BDF handles this efficiently; Adams needs many more steps.
tspan4 : np_linspace(0.0, 0.5, 200)$
/* BDF (stiff solver) */
res_bdf : np_cvode([-1000*y + 3000 - 2000*exp(-t)],
[t, y], [0.0], tspan4, 1e-10, 1e-10, bdf)$
/* Adams (non-stiff) for comparison */
res_adams : np_cvode([-1000*y + 3000 - 2000*exp(-t)],
[t, y], [0.0], tspan4, 1e-10, 1e-10, adams)$
print("Both methods converge to the same solution")$
t4 : np_slice(res_bdf, all, 0)$
y_bdf : np_slice(res_bdf, all, 1)$
y_adams : np_slice(res_adams, all, 1)$
print("Max difference:", np_max(np_abs(np_sub(y_bdf, y_adams))))$
Both methods converge to the same solution Max difference: 5.690126148039099e-10
/* Exact solution */
y_exact4(t) := 3 - 0.998*exp(-1000*t) - 2.002*exp(-t)$
ax_draw2d(
color="#2563eb", line_width=2, name="BDF",
lines(t4, y_bdf),
color="#dc2626", dash="dot", name="Exact",
explicit(y_exact4(t), t, 0, 0.5),
title="Stiff ODE: BDF Method", xlabel="t", ylabel="y",
grid=true, showlegend=true
)$
5. Coupled System — Brusselator¶
A classic chemical oscillator: $$\dot{x} = A + x^2 y - (B+1)x, \quad \dot{y} = Bx - x^2 y$$
With $A = 1$, $B = 3$, the system exhibits limit-cycle oscillations.
tspan5 : np_linspace(0.0, 40.0, 2000)$
A : 1.0$ B : 3.0$
res5 : np_cvode([A + x^2*y - (B+1)*x, B*x - x^2*y],
[t, x, y], [1.5, 3.0], tspan5)$
t5 : np_slice(res5, all, 0)$
x5 : np_slice(res5, all, 1)$
y5 : np_slice(res5, all, 2)$
ax_draw2d(
color="#2563eb", line_width=1.5, name="x(t)",
lines(t5, x5),
color="#f59e0b", line_width=1.5, name="y(t)",
lines(t5, y5),
title="Brusselator Oscillations", xlabel="t", ylabel="",
grid=true, showlegend=true
)$
/* Phase portrait shows the limit cycle */
ax_draw2d(
color="#7c3aed", line_width=1.5,
lines(x5, y5),
title="Brusselator Limit Cycle", xlabel="x", ylabel="y",
grid=true
)$
