Leveraging CVODE With the Robertson Problem
TL;DR. The three‑species Robertson kinetics form a classic stiff ODE system—perfect to demonstrate CVODE’s BDF/Newton toolbox and tolerance strategy.
The model
\[\begin{aligned} y_1′ & = -k_1 y_1 + k_3 y_2 y_3,\\ y_2′ & = k_1 y_1 – k_2 y_2^2 – k_3 y_2 y_3,\\ y_3′ & = k_2 y_2^2\,,\quad (k_1{=}0.04,\ k_2{=}10^4,\ k_3{=}3\cdot10^7) \end{aligned}\]Why stiffness matters
- Explicit methods choke on fast transients → tiny steps.
- Implicit BDF tracks slow manifolds efficiently.
CVODE setup (dense)
// Outline: init, tolerances, linear solver, integrate
SUNContext ctx; SUNContext_Create(0, &ctx);
N_Vector y = N_VNew_Serial(3, ctx);
NV_Ith_S(y,0)=1.0; NV_Ith_S(y,1)=0.0; NV_Ith_S(y,2)=0.0;
void *cvode = CVodeCreate(CV_BDF, ctx);
CVodeInit(cvode, f, 0.0, y);
CVodeSStolerances(cvode, 1e-6, 1e-10);
SUNMatrix A = SUNDenseMatrix(3,3,ctx);
SUNLinearSolver LS = SUNLinSol_Dense(y, A, ctx);
CVodeSetLinearSolver(cvode, LS, A);
// ... integrate, query stats ...Common pitfalls
Missing rpath → add -Wl,-rpath,$HOME/sundials-install/lib. Context required (v6+) → use SUNContext_Create and pass to vectors/solvers.