Stiff ODE Example — Robertson Chemical Kinetics¶
This example solves Robertson’s benchmark stiff 3-species chemical kinetics
system using differentialequations.ode_stiff() (Radau implicit Runge-Kutta)
and differentialequations.ivpgb() (Gear’s BDF method):
\[\begin{split}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\times10^7\,y_2^2 \\
y_3' &= 3\times10^7\,y_2^2\end{split}\]
with \(y_0 = [1, 0, 0]\) over \(t \in [0, 10^{11}]\).
Example Code¶
"""IMSL stiff ODE example: Robertson's chemical kinetics problem.
Demonstrates the use of ode_stiff() (Radau implicit RK) and ivpgb() (BDF)
for Robertson's stiff 3-species chemical kinetics system:
y1' = -0.04*y1 + 1e4*y2*y3
y2' = 0.04*y1 - 1e4*y2*y3 - 3e7*y2^2
y3' = 3e7*y2^2
with initial conditions y0 = [1, 0, 0] integrated over t in [0, 1e11].
Outputs:
- Table printed to stdout
- SVG plot saved to test_output/example_imsl_stiff_ode.svg
"""
from __future__ import annotations
from pathlib import Path
from typing import Dict
import matplotlib.pyplot as plt
import numpy as np
from differentialequations import ode_stiff, ivpgb
def robertson(t: float, y: np.ndarray) -> list:
"""Robertson stiff chemical kinetics RHS.
Args:
t (float): Time (unused, autonomous system).
y (np.ndarray): State vector [y1, y2, y3].
Returns:
list: Derivatives [y1', y2', y3'].
"""
y1, y2, y3 = y
dy1 = -0.04 * y1 + 1.0e4 * y2 * y3
dy2 = 0.04 * y1 - 1.0e4 * y2 * y3 - 3.0e7 * y2 ** 2
dy3 = 3.0e7 * y2 ** 2
return [dy1, dy2, dy3]
def run_demo_imsl_stiff_ode() -> Dict[str, object]:
"""Run IMSL stiff ODE example: Robertson chemical kinetics.
Solves Robertson's benchmark stiff ODE using both ode_stiff (Radau) and
ivpgb (BDF) and compares final-state results on a log time scale.
Args:
None
Returns:
Dict[str, object]: Result dict with keys ``t`` (ndarray),
``y`` (ndarray), and ``plot_path`` (str).
"""
y0 = [1.0, 0.0, 0.0]
t_span = (0.0, 1.0e11)
result_radau = ode_stiff(robertson, t_span, y0, rtol=1e-8, atol=1e-10)
result_bdf = ivpgb(robertson, t_span, y0, rtol=1e-8, atol=1e-10)
r = result_radau
yf = r.y[:, -1]
print("\nIMSL Stiff ODE Example: Robertson Chemical Kinetics")
print("=" * 60)
print(f"{'Parameter':<30} {'Radau (ode_stiff)':>20} {'BDF (ivpgb)':>12}")
print("-" * 60)
print(f"{'t_final':<30} {r.t[-1]:>20.4e} {result_bdf.t[-1]:>12.4e}")
print(f"{'y1 final':<30} {yf[0]:>20.8f} {result_bdf.y[0,-1]:>12.8f}")
print(f"{'y2 final (x 1e4)':<30} {yf[1]*1e4:>20.6e} {result_bdf.y[1,-1]*1e4:>12.6e}")
print(f"{'y3 final':<30} {yf[2]:>20.8f} {result_bdf.y[2,-1]:>12.8f}")
print(f"{'Steps taken':<30} {r.n_steps:>20d} {result_bdf.n_steps:>12d}")
print(f"{'Mass conservation y1+y3':<30} {yf[0]+yf[2]:>20.8f}")
print("=" * 60)
output_dir = Path("test_output")
output_dir.mkdir(parents=True, exist_ok=True)
plot_path = output_dir / "example_imsl_stiff_ode.svg"
t = r.t
y = r.y
t_safe = np.where(t > 0, t, np.finfo(float).tiny)
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
axes[0].semilogx(t_safe, y[0], color="#0369a1", linewidth=2.0, label=r"$y_1$")
axes[0].semilogx(t_safe, y[1] * 1e4, color="#e11d48", linewidth=1.8,
linestyle="--", label=r"$y_2 \times 10^4$")
axes[0].semilogx(t_safe, y[2], color="#16a34a", linewidth=1.8,
linestyle=":", label=r"$y_3$")
axes[0].set_xlabel("Time t")
axes[0].set_ylabel("Species concentration")
axes[0].set_title("Robertson Kinetics — Radau (ode_stiff)")
axes[0].legend()
axes[0].grid(True, which="both", alpha=0.3)
t2 = result_bdf.t
y2 = result_bdf.y
t2_safe = np.where(t2 > 0, t2, np.finfo(float).tiny)
axes[1].semilogx(t2_safe, y2[0], color="#0369a1", linewidth=2.0, label=r"$y_1$")
axes[1].semilogx(t2_safe, y2[1] * 1e4, color="#e11d48", linewidth=1.8,
linestyle="--", label=r"$y_2 \times 10^4$")
axes[1].semilogx(t2_safe, y2[2], color="#16a34a", linewidth=1.8,
linestyle=":", label=r"$y_3$")
axes[1].set_xlabel("Time t")
axes[1].set_ylabel("Species concentration")
axes[1].set_title("Robertson Kinetics — BDF (ivpgb)")
axes[1].legend()
axes[1].grid(True, which="both", alpha=0.3)
fig.suptitle("IMSL Stiff ODE: Robertson Chemical Kinetics", fontsize=13, fontweight="bold")
fig.tight_layout()
fig.savefig(plot_path, format="svg")
plt.close(fig)
return {"t": t, "y": y, "plot_path": str(plot_path)}
if __name__ == "__main__":
run_demo_imsl_stiff_ode()
Plot Output¶
Species concentrations versus log time. Left panel uses ode_stiff
(Radau); right panel uses ivpgb (BDF). Note that \(y_2\) is
scaled by \(10^4\) for visibility.¶
Console Output¶
IMSL Stiff ODE Example: Robertson Chemical Kinetics
============================================================
Parameter Radau (ode_stiff) BDF (ivpgb)
------------------------------------------------------------
t_final 1.0000e+11 1.0000e+11
y1 final 0.00000002 0.00000002
y2 final (x 1e4) 8.333359e-10 8.405940e-10
y3 final 0.99999998 0.99999998
Steps taken 787 979
Mass conservation y1+y3 1.00000000
============================================================