Runnable example

Robertson’s problem, solved three ways.

Robertson’s autocatalytic reaction is a standard test for stiffness. Its rate constants span nine orders of magnitude, so the fast transient is over in microseconds while the slow part runs out to $10^{7}$ seconds. Below, the same system is solved three ways with Reineur Solver, a numerical library for C#.

\begin{aligned} \dot{y}_1 &= -0.04\,y_1 + 10^{4}\,y_2 y_3 \\ \dot{y}_2 &= 0.04\,y_1 - 10^{4}\,y_2 y_3 - 3\times10^{7} y_2^{2} \\ \dot{y}_3 &= 3\times10^{7} y_2^{2} \end{aligned}

In the differential-algebraic form the third equation is replaced by the conservation constraint $y_1 + y_2 + y_3 = 1$ , which makes the mass matrix singular.

The code is editable. Replace it with any problem and press Run. It compiles and runs on a server, so there is nothing to install and no account to create.

Mass-matrix form. The third equation is the algebraic constraint y₁ + y₂ + y₃ = 1, so M is singular and the system is an index-1 DAE.

C#, editable16 lines

double[] robertson_f(double t, double[] y) =>
    [(-0.04 * y[0] + 1e4 * y[1] * y[2]),
     (0.04 * y[0] - 1e4 * y[1] * y[2] - 3e7 * y[1]*y[1]),
     y[0] + y[1] + y[2] - 1.0];
double[,] mass_f = Diag([1, 1, 0]);
double[] y0 = [1.0, 0.0, 0.0];
(ColVec T, Matrix Y) = Ode43a(robertson_f, mass_f, y0, [0, 1e7]);
// Plot the result
Y[.., 1] = 1e4*Y[.., 1];
SemiLogx(T, Y);
Xlabel("Time t"); Ylabel("Soluton y");
Legend(["y_1", "1e4*y_2", "y_3"], MiddleLeft);
Title("Solution of Robertson's ODE with ODE43a");
SaveAs("Robertson-ODE-given-points-Ode43a.png");

This chapter in the book

Result

Press Run to solve this system.
Output and any plot will appear here.

Method vs method

Comparing methods on the same problem

Most results in numerical methods come from comparing two approaches on a problem chosen to show the difference between them. Both sides below actually run. Edit either one and press Run.

Why can't you just use a general-purpose solver?

Van der Pol at mu = 100 is stiff: the solution has a fast transient riding on a slow limit cycle. Both solvers below get the right answer. Watch what it costs the explicit one.

Explicit Runge-Kutta

Ode45

double[] vdp(double t, double[] y) =>
    [y[1], 100 * (1 - y[0]*y[0]) * y[1] - y[0]];
var opts = Odeset(Stats: true);
(ColVec T, Matrix Y) = Ode45(vdp, [2, 0], [0, 300], opts);
Plot(T, Y[.., 0], "b");
Xlabel("t"); Ylabel("y_1");
Title("Van der Pol, mu = 100 - Ode45");
SaveAs("vdp_Ode45.png");

Not run yet

Stiff solver

Ode43s

double[] vdp(double t, double[] y) =>
    [y[1], 100 * (1 - y[0]*y[0]) * y[1] - y[0]];
var opts = Odeset(Stats: true);
(ColVec T, Matrix Y) = Ode43s(vdp, [2, 0], [0, 300], opts);
Plot(T, Y[.., 0], "b");
Xlabel("t"); Ylabel("y_1");
Title("Van der Pol, mu = 100 - Ode43s");
SaveAs("vdp_Ode43s.png");

Not run yet