#equation-solver #differential #numerical-methods #solver #numerical

bin+lib russell_ode

Solvers for ordinary differential equations and differential algebraic equations

17 stable releases

1.9.0 Sep 28, 2024
1.8.0 Sep 25, 2024
1.6.2 Jul 23, 2024
1.6.0 May 21, 2024
0.9.2 Mar 17, 2024

#69 in Math

MIT license

3MB
41K SLoC

Russell ODE - Solvers for ordinary differential equations and differential algebraic equations

documentation

This crate is part of Russell - Rust Scientific Library

Contents

Introduction

This library implements (in native Rust) solvers for ordinary differential equations (ODEs) and differential algebraic equations (DAEs). Specifically, it implements several explicit Runge-Kutta methods (e.g., Dormand-Prince formulae) and two implicit Runge-Kutta methods, namely the Backward Euler and the Radau IIA of fifth-order (aka Radau5). The Radau5 solver can solve DAEs of Index-1 by accepting the so-called mass matrix.

The code in this library is based on Hairer-NΓΈrsett-Wanner books and respective Fortran codes (see references [1] and [2]). However, the implementations of Dormand-Prince 5(4) and Dormand-Prince 8(5,3) are different from the Fortran counterparts. The code for Radau5 follows closely reference [2]; however, some minor differences are considered. Despite the coding differences, the numeric results match the Fortran results quite well.

The recommended methods are:

  • DoPri5 for ODE systems and non-stiff problems using moderate tolerances
  • DoPri8 for ODE systems and non-stiff problems using strict tolerances
  • Radau5 for ODE and DAE systems, possibly stiff, with moderate to strict tolerances

The ODE/DAE system can be easily defined using the System data structure; see the examples below.

Documentation

References

  1. Hairer E, NΓΈrsett, SP, Wanner G (2008) Solving Ordinary Differential Equations I. Non-stiff Problems. Second Revised Edition. Corrected 3rd printing 2008. Springer Series in Computational Mathematics, 528p
  2. Hairer E, Wanner G (2002) Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems. Second Revised Edition. Corrected 2nd printing 2002. Springer Series in Computational Mathematics, 614p
  3. Kreyszig, E (2011) Advanced engineering mathematics; in collaboration with Kreyszig H, Edward JN 10th ed 2011, Hoboken, New Jersey, Wiley

Installation

This crate depends on some non-rust high-performance libraries. See the main README file for the steps to install these dependencies.

Setting Cargo.toml

Crates.io

πŸ‘† Check the crate version and update your Cargo.toml accordingly:

[dependencies]
russell_ode = "*"

Optional features

The following (Rust) features are available:

  • intel_mkl: Use Intel MKL instead of OpenBLAS
  • local_suitesparse: Use a locally compiled version of SuiteSparse
  • with_mumps: Enable the MUMPS solver (locally compiled)

Note that the main README file presents the steps to compile the required libraries according to each feature.

🌟 Examples

This section illustrates how to use russell_ode. See also:

Simple ODE with a single equation

Solve the simple ODE with Dormand-Prince 8(5,3):

      dy
y' = β€”β€” = 1   with   y(x=0)=0    thus   y(x) = x
      dx

See the code simple_ode_single_equation.rs; reproduced below:

use russell_lab::{vec_approx_eq, StrError, Vector};
use russell_ode::prelude::*;

fn main() -> Result<(), StrError> {
    // system
    let ndim = 1;
    let system = System::new(ndim, |f, x, y, _args: &mut NoArgs| {
        f[0] = x + y[0];
        Ok(())
    });

    // solver
    let params = Params::new(Method::DoPri8);
    let mut solver = OdeSolver::new(params, system)?;

    // initial values
    let x = 0.0;
    let mut y = Vector::from(&[0.0]);

    // solve from x = 0 to x = 1
    let x1 = 1.0;
    let mut args = 0;
    solver.solve(&mut y, x, x1, None, &mut args)?;
    println!("y =\n{}", y);

    // check the results
    let y_ana = Vector::from(&[f64::exp(x1) - x1 - 1.0]);
    vec_approx_eq(&y, &y_ana, 1e-7);

    // print stats
    println!("{}", solver.stats());
    Ok(())
}

The output looks like:

y =
β”Œ                   ┐
β”‚ 0.718281815054018 β”‚
β””                   β”˜
DoPri8: Dormand-Prince method (explicit, order 8(5,3), embedded)
Number of function evaluations   = 84
Number of performed steps        = 7
Number of accepted steps         = 7
Number of rejected steps         = 0
Last accepted/suggested stepsize = 0.40139999999999776
Max time spent on a step         = 7.483Β΅s
Total time                       = 65.462Β΅s

Simple system with mass matrix

Solve with Radau5:

y0' + y1'     = -y0 + y1
y0' - y1'     =  y0 + y1
          y2' = 1/(1 + x)

y0(0) = 1,  y1(0) = 0,  y2(0) = 0

Thus:

M y' = f(x, y)

with:

    β”Œ          ┐       β”Œ           ┐
    β”‚  1  1  0 β”‚       β”‚ -y0 + y1  β”‚
M = β”‚  1 -1  0 β”‚   f = β”‚  y0 + y1  β”‚
    β”‚  0  0  1 β”‚       β”‚ 1/(1 + x) β”‚
    β””          β”˜       β””           β”˜

The Jacobian matrix is:

         β”Œ          ┐
    df   β”‚ -1  1  0 β”‚
J = β€”β€” = β”‚  1  1  0 β”‚
    dy   β”‚  0  0  0 β”‚
         β””          β”˜

The analytical solution is:

y0(x) = cos(x)
y1(x) = -sin(x)
y2(x) = log(1 + x)

Reference: Numerical Solution of Differential-Algebraic Equations: Solving Systems with a Mass Matrix.

See the code simple_system_with_mass.rs; reproduced below:

use russell_lab::{vec_approx_eq, StrError, Vector};
use russell_ode::prelude::*;
use russell_sparse::{CooMatrix, Sym};

fn main() -> Result<(), StrError> {
    // ODE system
    let ndim = 3;
    let jac_nnz = 4;
    let mut system = System::new(ndim, |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| {
        f[0] = -y[0] + y[1];
        f[1] = y[0] + y[1];
        f[2] = 1.0 / (1.0 + x);
        Ok(())
    });


    // function to compute the Jacobian matrix
    let symmetric = Sym::No;
    system.set_jacobian(
        Some(jac_nnz),
        symmetric,
        move |jj: &mut CooMatrix, alpha: f64, _x: f64, _y: &Vector, _args: &mut NoArgs| {
            jj.reset();
            jj.put(0, 0, alpha * (-1.0))?;
            jj.put(0, 1, alpha * (1.0))?;
            jj.put(1, 0, alpha * (1.0))?;
            jj.put(1, 1, alpha * (1.0))?;
            Ok(())
        },
    )?;

    // mass matrix
    let mass_nnz = 5;
    system.set_mass(Some(mass_nnz), symmetric, |mm: &mut CooMatrix| {
        mm.put(0, 0, 1.0).unwrap();
        mm.put(0, 1, 1.0).unwrap();
        mm.put(1, 0, 1.0).unwrap();
        mm.put(1, 1, -1.0).unwrap();
        mm.put(2, 2, 1.0).unwrap();
    })?;

    // solver
    let params = Params::new(Method::Radau5);
    let mut solver = OdeSolver::new(params, system)?;

    // initial values
    let x = 0.0;
    let mut y = Vector::from(&[1.0, 0.0, 0.0]);

    // solve from x = 0 to x = 20
    let x1 = 20.0;
    let mut args = 0;
    solver.solve(&mut y, x, x1, None, &mut args)?;
    println!("y =\n{}", y);

    // check the results
    let y_ana = Vector::from(&[f64::cos(x1), -f64::sin(x1), f64::ln(1.0 + x1)]);
    vec_approx_eq(&y, &y_ana, 1e-3);

    // print stats
    println!("{}", solver.stats());
    Ok(())
}

The output looks like:

y =
β”Œ                     ┐
β”‚ 0.40864108577398206 β”‚
β”‚ -0.9136567566808349 β”‚
β”‚   3.044521229909652 β”‚
β””                     β”˜
Radau5: Radau method (Radau IIA) (implicit, order 5, embedded)
Number of function evaluations   = 203
Number of Jacobian evaluations   = 1
Number of factorizations         = 8
Number of lin sys solutions      = 52
Number of performed steps        = 46
Number of accepted steps         = 46
Number of rejected steps         = 0
Number of iterations (maximum)   = 2
Number of iterations (last step) = 1
Last accepted/suggested stepsize = 0.5055117216674699
Max time spent on a step         = 40.817Β΅s
Max time spent on the Jacobian   = 558ns
Max time spent on factorization  = 199.895Β΅s
Max time spent on lin solution   = 53.101Β΅s
Total time                       = 2.653323ms

Brusselator ODE

This example corresponds to Fig 16.4 on page 116 of Reference #1. The problem is defined in Eq (16.12) on page 116 of Reference #1.

The system is:

y0' = 1 - 4 y0 + y0Β² y1
y1' = 3 y0 - y0Β² y1

with  y0(x=0) = 3/2  and  y1(x=0) = 3

The Jacobian matrix is:

         β”Œ                     ┐
    df   β”‚ -4 + 2 y0 y1    y0Β² β”‚
J = β€”β€” = β”‚                     β”‚
    dy   β”‚  3 - 2 y0 y1   -y0Β² β”‚
         β””                     β”˜

Solving with DoPri8 -- 8(5,3) -- dense output

This is a system of two ODEs, well explained in Reference #1. This problem is solved with the DoPri8 method (it has a hybrid error estimator of 5th and 3rd order; see Reference #1).

This example also shows how to enable the dense output.

See the code brusselator_ode_dopri8.rs; reproduced below (without the plotting commands):

The output looks like this:

y_russell     = [0.4986435155366857, 4.596782273713258]
y_mathematica = [0.49863707126834783, 4.596780349452011]
DoPri8: Dormand-Prince method (explicit, order 8(5,3), embedded)
Number of function evaluations   = 647
Number of performed steps        = 45
Number of accepted steps         = 38
Number of rejected steps         = 7
Last accepted/suggested stepsize = 2.1617616186304227
Max time spent on a step         = 47.643Β΅s
Total time                       = 898.347Β΅s

A plot of the (dense) solution is shown below:

Brusselator results: DoPri8

Variable step sizes

This example solves the Brusselator ODE with variable step sizes for different tolerances. In this example, tol = abs_tol = rel_tol. The global error is the difference between Russell's results and Mathematica's results obtained with high accuracy. The Mathematica code is:

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
sys = GetNDSolveProblem["BrusselatorODE"];
sol = NDSolve[sys, Method -> "StiffnessSwitching", WorkingPrecision -> 32];
ref = First[FinalSolutions[sys, sol]]

See the code brusselator_ode_var_step.rs

The results are:

━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
       tol =  1.00E-02  1.00E-04  1.00E-06  1.00E-08
      Method     Error     Error     Error     Error
────────────────────────────────────────────────────
      Radau5   1.9E-03   7.9E-06   1.3E-07   3.2E-09
     Merson4   3.8E-02   2.1E-04   9.9E-06   7.1E-08
      DoPri5   8.0E-03   1.7E-04   1.8E-06   2.0E-08
      DoPri8   1.4E-02   6.4E-06   2.7E-07   4.2E-09
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

And the convergence plot is:

Brusselator results: var step

Fixed step sizes

This example solves the Brusselator ODE with fixed step sizes and explicit Runge-Kutta methods. The global error is also the difference between Russell and Mathematica as in the previous section.

See the code brusselator_ode_fix_step.rs

The results are:

━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
        h = 2.00E-01 1.00E-01 5.00E-02 1.00E-02 1.00E-03
     Method    Error    Error    Error    Error    Error
────────────────────────────────────────────────────────
        Rk2  6.1E-03  4.1E-03  1.5E-03  7.4E-05  7.8E-07
        Rk3  6.1E-03  8.2E-04  1.5E-04  1.7E-06  1.8E-09
      Heun3  1.1E-02  1.3E-03  1.7E-04  1.5E-06  1.5E-09
        Rk4  5.9E-03  2.8E-04  1.7E-05  2.7E-08  2.7E-12
     Rk4alt  7.0E-03  2.4E-04  9.3E-06  1.1E-08  9.9E-13
    MdEuler  4.1E-02  7.2E-03  2.0E-03  9.0E-05  9.2E-07
    Merson4  5.7E-04  6.3E-06  7.9E-07  1.7E-09  1.6E-13
 Zonneveld4  5.9E-03  2.8E-04  1.7E-05  2.7E-08  2.7E-12
  Fehlberg4  9.4E-03  8.1E-05  1.4E-06  1.7E-09  1.7E-13
     DoPri5  1.9E-03  6.1E-05  7.3E-07  4.7E-11  5.7E-14
    Verner6  3.6E-03  3.8E-05  3.5E-07  1.4E-11  4.1E-14
  Fehlberg7  2.5E-05  9.9E-08  6.3E-10  8.9E-15  8.9E-15
     DoPri8  3.9E-06  5.3E-10  5.8E-12  1.7E-14  2.0E-14
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

And the convergence plot is:

Brusselator results: fix step

Brusselator PDE

This example corresponds to Fig 10.4(a,b) on pages 250 and 251 of Reference #1. The problem is defined in Eqs (10.10-10.14) on pages 248 and 249 of Reference #1.

If second_book is true, this example corresponds to Fig 10.7 on page 151 of Reference #2. Also, in this case, the problem is defined in Eqs (10.15-10.16) on pages 151 and 152 of Reference #2.

While in the first book the boundary conditions are Neumann-type, in the second book the boundary conditions are periodic. Also the initial values in the second book are different that those in the first book.

The model is given by:

βˆ‚u                         βŽ› βˆ‚Β²u   βˆ‚Β²u ⎞
β€”β€”β€” = 1 - 4.4 u + uΒ² v + Ξ± ⎜ β€”β€”β€” + β€”β€”β€” ⎟ + I(t,x,y)
βˆ‚t                         ⎝ βˆ‚xΒ²   βˆ‚yΒ² ⎠

βˆ‚v                         βŽ› βˆ‚Β²v   βˆ‚Β²v ⎞
β€”β€”β€” =     3.4 u - uΒ² v + Ξ± ⎜ β€”β€”β€” + β€”β€”β€” ⎟
βˆ‚t                         ⎝ βˆ‚xΒ²   βˆ‚yΒ² ⎠

with:  t β‰₯ 0,  0 ≀ x ≀ 1,  0 ≀ y ≀ 1

where I(t,x,y) is the inhomogeneity function (second book) given by:

            ⎧ 5  if (x-0.3)Β²+(y-0.6)Β² ≀ 0.1Β² and t β‰₯ 1.1
I(t,x,y) =  ⎨
            ⎩ 0  otherwise

The first book considers the following Neumann boundary conditions:

βˆ‚u          βˆ‚v     
β€”β€”β€” = 0     β€”β€”β€” = 0
 β†’           β†’
βˆ‚n          βˆ‚n     

and the following initial conditions (first book):

u(t=0,x,y) = 0.5 + y    v(t=0,x,y) = 1 + 5 x

The second book considers periodic boundary conditions on u. However, here we assume periodic on u and v:

u(t, 0, y) = u(t, 1, y)
u(t, x, 0) = u(t, x, 1)
v(t, 0, y) = v(t, 1, y)   ← Not in the book
v(t, x, 0) = v(t, x, 1)   ← Not in the book

The second book considers the following initial conditions:

u(0, x, y) = 22 y pow(1 - y, 1.5)
v(0, x, y) = 27 x pow(1 - x, 1.5)

The scalar fields u(x, y) and v(x, y) are mapped over a rectangular grid with their discrete counterparts represented by:

pα΅’β±Ό(t) := u(t, xα΅’, yβ±Ό)
qα΅’β±Ό(t) := v(t, xα΅’, yβ±Ό)

Thus ndim = 2 npointΒ² with npoint being the number of points along the x or y line.

The second partial derivatives over x and y (Laplacian) are approximated using the Finite Differences Method (FDM).

The pα΅’β±Ό and qα΅’β±Ό values are mapped onto the vectors U and V as follows:

pα΅’β±Ό β†’ Uβ‚˜
qα΅’β±Ό β†’ Vβ‚˜

with m = i + j nx

Then, they are stored in a single vector Y:

    β”Œ   ┐
    β”‚ U β”‚
Y = β”‚   β”‚
    β”‚ V β”‚
    β””   β”˜

Thus:

Uβ‚˜ = Yβ‚˜  and  Vβ‚˜ = Yβ‚›β‚Šβ‚˜

where  0 ≀ m ≀ s - 1  and (shift)  s = npointΒ²

In terms of components, we can write:

      ⎧ Uₐ    if a < s
Yₐ =  ⎨
      ⎩ Vₐ₋ₛ  if a β‰₯ s

where  0 ≀ a ≀ ndim - 1  and  ndim = 2 s

The components of the resulting system of equations are defined by: (the prime indicates time-derivative; no summation over repeated indices):

Uβ‚˜' = 1 - 4.4 Uβ‚˜ + Uβ‚˜Β² Vβ‚˜ + Ξ£ Aβ‚˜β‚– Uβ‚–
                           k
Vβ‚˜' =     3.4 Uβ‚˜ - Uβ‚˜Β² Vβ‚˜ + Ξ£ Aβ‚˜β‚– Uβ‚–
                           k

where Aβ‚˜β‚– are the elements of the discrete Laplacian matrix

The components to build the Jacobian matrix are: (no summation over repeated indices):

βˆ‚Uβ‚˜'
β€”β€”β€”β€” = -4.4 Ξ΄β‚˜β‚™ + 2 Uβ‚˜ Ξ΄β‚˜β‚™ Vβ‚˜ + Aβ‚˜β‚™
βˆ‚Uβ‚™

βˆ‚Uβ‚˜'
β€”β€”β€”β€” = Uβ‚˜Β² Ξ΄β‚˜β‚™
βˆ‚Vβ‚™

βˆ‚Vβ‚˜'
β€”β€”β€”β€” = 3.4 Ξ΄β‚˜β‚™ - 2 Uβ‚˜ Ξ΄β‚˜β‚™ Vβ‚˜
βˆ‚Uβ‚™

βˆ‚Vβ‚˜'
β€”β€”β€”β€” = -Uβ‚˜Β² Ξ΄β‚˜β‚™ + Aβ‚˜β‚™
βˆ‚Vβ‚™

where Ξ΄β‚˜β‚™ is the Kronecker delta

With Fₐ := βˆ‚Yₐ/βˆ‚t, the components of the Jacobian matrix can be "assembled" as follows:

      ⎧  ⎧ βˆ‚Uₐ'/βˆ‚Uβ‚‘      if e < s
      β”‚  ⎨                          if a < s
βˆ‚Fₐ   β”‚  ⎩ βˆ‚Uₐ'/βˆ‚Vβ‚‘β‚‹β‚›    if e β‰₯ s
β€”β€”β€” = ⎨
βˆ‚Yβ‚‘   β”‚  ⎧ βˆ‚Vₐ₋ₛ'/βˆ‚Uβ‚‘    if e < s
      β”‚  ⎨                          if a β‰₯ s
      ⎩  ⎩ βˆ‚Vₐ₋ₛ'/βˆ‚Vβ‚‘β‚‹β‚›  if e β‰₯ s

where  0 ≀ a ≀ ndim - 1  and  0 ≀ e ≀ ndim - 1

First book

We solve this problem with Radau5. The approximated solution is generated with npoint = 21 and is implemented in brusselator_pde_radau5.rs. This code will generate a series of files, one for each (dense) output time with h_out = 0.5.

The results can then be plotted with brusselator_pde_plot.rs

The results are shown below for the U field:

brusselator_pde_radau5_u.svg

And below for the V field:

brusselator_pde_radau5_v.svg

These figures compare well with the corresponding ones on pages 250 and 251 of Reference #1.

The computations with russell are also compared with values obtained with Mathematica. The verification is implemented in test_radau5_brusselator_pde.

The figure below shows the russell (black dashed lines) and Mathematica (red solid lines) results for the U field:

test_radau5_brusselator_pde_u.svg

The figure below shows the russell (black dashed lines) and Mathematica (red solid lines) results for the V field:

test_radau5_brusselator_pde_v.svg

Second book

For the problem in the second book, we run brusselator_pde_radau5_2nd.rs with npoint = 129 and h_out = 1.0.

The results are shown below for the U field:

brusselator_pde_radau5_2nd_u.jpg

And below for the V field:

brusselator_pde_radau5_2nd_v.jpg

The code brusselator_pde_2nd_comparison.rs compares russell results with Mathematica results.

The figure below shows the russell (black dashed lines) and Mathematica (red solid lines) results for the U field:

comparison U

The figure below shows the russell (black dashed lines) and Mathematica (red solid lines) results for the V field:

comparison V

Arenstorf orbits

This example corresponds to Fig 0.1 on page 130 of Reference #1. The problem is defined in Eqs (0.1) and (0.2) on page 129 and 130 of Reference #1.

From Hairer-NΓΈrsett-Wanner:

"(...) an example from Astronomy, the restricted three body problem. (...) two bodies of masses ΞΌ' = 1 βˆ’ ΞΌ and ΞΌ in circular rotation in a plane and a third body of negligible mass moving around in the same plane. (...)"

The system equations are:

y0'' = y0 + 2 y1' - ΞΌ' (y0 + ΞΌ) / d0 - ΞΌ (y0 - ΞΌ') / d1
y1'' = y1 - 2 y0' - ΞΌ' y1 / d0 - ΞΌ y1 / d1

With the assignments:

y2 := y0'  β‡’  y2' = y0''
y3 := y1'  β‡’  y3' = y1''

We obtain a 4-dim problem:

f0 := y0' = y2
f1 := y1' = y3
f2 := y2' = y0 + 2 y3 - ΞΌ' (y0 + ΞΌ) / d0 - ΞΌ (y0 - ΞΌ') / d1
f3 := y3' = y1 - 2 y2 - ΞΌ' y1 / d0 - ΞΌ y1 / d1

See the code arenstorf_dopri8.rs

The code output is:

y_russell     = [0.9943002573065823, 0.000505756322923528, 0.07893182893575335, -1.9520617089599261]
y_mathematica = [0.9939999999999928, 2.4228439406717e-14, 3.6631563591513e-12, -2.0015851063802006]
DoPri8: Dormand-Prince method (explicit, order 8(5,3), embedded)
Number of function evaluations   = 870
Number of performed steps        = 62
Number of accepted steps         = 47
Number of rejected steps         = 15
Last accepted/suggested stepsize = 0.005134142939114363
Max time spent on a step         = 10.538Β΅s
Total time                       = 1.399021ms

The results are plotted below:

Arenstorf Orbit: DoPri8

Hairer-Wanner Equation (1.1)

This example corresponds to Fig 1.1 and Fig 1.2 on page 2 of Reference #2. The problem is defined in Eq (1.1) on page 2 of Reference #2.

The system is:

y0' = -50 (y0 - cos(x))

with  y0(x=0) = 0

The Jacobian matrix is:

    df   β”Œ     ┐
J = β€”β€” = β”‚ -50 β”‚
    dy   β””     β”˜

This example illustrates the instability of the forward Euler method with step sizes above the stability limit. The equation is (reference # 2, page 2):

This example also shows how to enable the output of accepted steps.

See the code hairer_wanner_eq1.rs

The results are show below:

Hairer-Wanner Eq(1.1)

Robertson's Equation

This example corresponds to Fig 1.3 on page 4 of Reference #2. The problem is defined in Eq (1.4) on page 3 of Reference #2.

The system is:

y0' = -0.04 y0 + 1.0e4 y1 y2
y1' =  0.04 y0 - 1.0e4 y1 y2 - 3.0e7 y1Β²
y2' =                          3.0e7 y1Β²

with  y0(0) = 1, y1(0) = 0, y2(0) = 0

This example illustrates the Robertson's equation. In this problem DoPri5 uses many steps (about 200). On the other hand, Radau5 solves the problem with 17 accepted steps.

This example also shows how to enable the output of accepted steps.

See the code robertson.rs

The solution is approximated with Radau5 and DoPri5 using two sets of tolerances.

The output is:

Radau5: Radau method (Radau IIA) (implicit, order 5, embedded)
Number of function evaluations   = 88
Number of Jacobian evaluations   = 8
Number of factorizations         = 15
Number of lin sys solutions      = 24
Number of performed steps        = 17
Number of accepted steps         = 15
Number of rejected steps         = 1
Number of iterations (maximum)   = 2
Number of iterations (last step) = 1
Last accepted/suggested stepsize = 0.8160578540023586
Max time spent on a step         = 117.916Β΅s
Max time spent on the Jacobian   = 1.228Β΅s
Max time spent on factorization  = 199.151Β΅s
Max time spent on lin solution   = 56.767Β΅s
Total time                       = 1.922108ms

Tol = 1e-2
DoPri5: Dormand-Prince method (explicit, order 5(4), embedded)
Number of function evaluations   = 1585
Number of performed steps        = 264
Number of accepted steps         = 209
Number of rejected steps         = 55
Last accepted/suggested stepsize = 0.0017137362591141277
Max time spent on a step         = 2.535Β΅s
Total time                       = 2.997516ms

Tol = 1e-3
DoPri5: Dormand-Prince method (explicit, order 5(4), embedded)
Number of function evaluations   = 1495
Number of performed steps        = 249
Number of accepted steps         = 205
Number of rejected steps         = 44
Last accepted/suggested stepsize = 0.0018175194753331549
Max time spent on a step         = 1.636Β΅s
Total time                       = 3.705391ms

The results are illustrated below:

Robertson's Equation - Solution

The step sizes from the DoPri solution with Tol = 1e-2 are illustrated below:

Robertson's Equation - Step Sizes

Van der Pol's Equation

This example corresponds Eq (1.5') on page 5 of Reference #2.

The system is:

y0' = y1
y1' = ((1.0 - y[0] * y[0]) * y[1] - y[0]) / Ξ΅

where Ξ΅ defines the stiffness of the problem + conditions (equation + initial conditions + step size + method).

DoPri5

This example corresponds to Fig 2.6 on page 23 of Reference #2. The problem is defined in Eq (1.5') on page 5 of Reference #2.

This example illustrated the stiffness of the Van der Pol problem with Ξ΅ = 0.003. In this example, DoPri5 with Tol = 1e-3 is used.

This example also shows how to enable the stiffness detection.

See the code van_der_pol_dopri5.rs

The output is given below:

y =
β”Œ                     ┐
β”‚   1.819918013289893 β”‚
β”‚ -0.7863062155442466 β”‚
β””                     β”˜
DoPri5: Dormand-Prince method (explicit, order 5(4), embedded)
Number of function evaluations   = 3133
Number of performed steps        = 522
Number of accepted steps         = 498
Number of rejected steps         = 24
Last accepted/suggested stepsize = 0.004363549192919735
Max time spent on a step         = 2.558Β΅s
Total time                       = 1.715917ms

The results are show below:

Van der Pol's Equation - DoPri5

The figure's red dashed lines mark the moment when stiffness has been detected first. The stiffness is confirmed after 15 accepted steps with repeated stiffness thresholds being reached. The positive thresholds are counted when h·ρ becomes greater than the corresponding factor·max(h·ρ)---the value on the stability limit (3.3 for DoPri5; factor ~= 0.976). Note that ρ is the approximation of the dominant eigenvalue of the Jacobian. After 6 accepted steps, if the thresholds are not reached, the stiffness detection flag is set to false.

Radau5

This example corresponds to Fig 8.1 on page 125 of Reference #2. The problem is defined in Eq (1.5') on page 5 of Reference #2.

This example uses a smaller Ξ΅ = 1e-6, making the problem + conditions much more stiff. It is solved with the Radau5 solver, which can handle stiff problems quite well. Note that DoPri5 would not solve this problem with such small Ξ΅, unless a very high number of steps (and oder configurations) were considered.

The output is given below:

y =
β”Œ                    ┐
β”‚ 1.7061626037853908 β”‚
β”‚ -0.892799551109113 β”‚
β””                    β”˜
Radau5: Radau method (Radau IIA) (implicit, order 5, embedded)
Number of function evaluations   = 2237
Number of Jacobian evaluations   = 160
Number of factorizations         = 252
Number of lin sys solutions      = 663
Number of performed steps        = 280
Number of accepted steps         = 241
Number of rejected steps         = 7
Number of iterations (maximum)   = 6
Number of iterations (last step) = 3
Last accepted/suggested stepsize = 0.2466642610579514
Max time spent on a step         = 136.487Β΅s
Max time spent on the Jacobian   = 949ns
Max time spent on factorization  = 223.917Β΅s
Max time spent on lin solution   = 4.010536ms
Total time                       = 37.8697ms

The results are show below:

Van der Pol's Equation - Radau5

One-transistor amplifier

This example corresponds to Fig 1.3 on page 377 and Fig 1.4 on page 379 of Reference #2. The problem is defined in Eq (1.14) on page 377 of Reference #2.

This is a differential-algebraic problem modelling the nodal voltages of a one-transistor amplifier.

The DAE is expressed in the so-called mass-matrix form (ndim = 5):

M y' = f(x, y)

with: y0(0)=0, y1(0)=Ub/2, y2(0)=Ub/2, y3(0)=Ub, y4(0)=0

where the elements of the right-hand side function are:

f0 = (y0 - ue) / R
f1 = (2 y1 - UB) / S + Ξ³ g12
f2 = y2 / S - g12
f3 = (y3 - UB) / S + Ξ± g12
f4 = y4 / S

with:

ue = A sin(Ο‰ x)
g12 = Ξ² (exp((y1 - y2) / UF) - 1)

Compared to Eq (1.14), we set all resistances Rα΅’ to S, except the first one (R := Rβ‚€).

The mass matrix is:

    β”Œ                     ┐
    β”‚ -C1  C1             β”‚
    β”‚  C1 -C1             β”‚
M = β”‚         -C2         β”‚
    β”‚             -C3  C3 β”‚
    β”‚              C3 -C3 β”‚
    β””                     β”˜

and the Jacobian matrix is:

    β”Œ                                           ┐
    β”‚ 1/R                                       β”‚
    β”‚       2/S + Ξ³ h12      -Ξ³ h12             β”‚
J = β”‚              -h12   1/S + h12             β”‚
    β”‚             Ξ± h12      -Ξ± h12             β”‚
    β”‚                                 1/S       β”‚
    β”‚                                       1/S β”‚
    β””                                           β”˜

with:

h12 = Ξ² exp((y1 - y2) / UF) / UF

Note: In this library, only Radau5 can solve such DAE.

See the code amplifier1t_radau5.rs

The output is given below:

y_russell     = [-0.022267, 3.068709, 2.898349, 1.499405, -1.735090]
y_mathematica = [-0.022267, 3.068709, 2.898349, 1.499439, -1.735057]
Radau5: Radau method (Radau IIA) (implicit, order 5, embedded)
Number of function evaluations   = 6007
Number of Jacobian evaluations   = 480
Number of factorizations         = 666
Number of lin sys solutions      = 1840
Number of performed steps        = 666
Number of accepted steps         = 481
Number of rejected steps         = 39
Number of iterations (maximum)   = 6
Number of iterations (last step) = 1
Last accepted/suggested stepsize = 0.00007705697843645314
Max time spent on a step         = 55.281Β΅s
Max time spent on the Jacobian   = 729ns
Max time spent on factorization  = 249.11Β΅s
Max time spent on lin solution   = 241.201Β΅s
Total time                       = 97.951021ms

The results are plotted below:

One-transistor Amplifier - Radau5

PDE: discrete Laplacian operator in 2D

For convenience (e.g., in benchmarks), russell_ode implements a discrete Laplacian operator (2D) based on the Finite Differences Method.

This operator can be used to solve simple partial differential equation (PDE) problems.

Laplace equation

Approximate (with the Finite Differences Method, FDM) the solution of

βˆ‚Β²Ο•     βˆ‚Β²Ο•
β€”β€”β€”  +  β€”β€”β€” = 0
βˆ‚xΒ²     βˆ‚yΒ²

on a (1.0 Γ— 1.0) rectangle with the following essential (Dirichlet) boundary conditions:

left:    Ο•(0.0, y) = 50.0
right:   Ο•(1.0, y) =  0.0
bottom:  Ο•(x, 0.0) =  0.0
top:     Ο•(x, 1.0) = 50.0

See the code pde_laplace_equation.rs

The results are illustrated below:

Laplace equation

Poisson equation 1

Approximate (with the Finite Differences Method, FDM) the solution of

βˆ‚Β²Ο•   βˆ‚Β²Ο•
β€”β€”β€” + β€”β€”β€” = 2 x (y - 1) (y - 2 x + x y + 2) exp(x - y)
βˆ‚xΒ²   βˆ‚yΒ²

on a (1.0 Γ— 1.0) square with the homogeneous boundary conditions.

The analytical solution is:

Ο•(x, y) = x y (x - 1) (y - 1) exp(x - y)

See the code test_pde_poisson_1.rs

The results are illustrated below:

Poisson equation 1

Poisson equation 2

Approximate (with the Finite Differences Method, FDM) the solution of

βˆ‚Β²Ο•   βˆ‚Β²Ο•
β€”β€”β€” + β€”β€”β€” = - π² y sin(Ο€ x)
βˆ‚xΒ²   βˆ‚yΒ²

on a (1.0 Γ— 1.0) square with the following essential boundary conditions:

left:    Ο•(0.0, y) = 0.0
right:   Ο•(1.0, y) = 0.0
bottom:  Ο•(x, 0.0) = 0.0
top:     Ο•(x, 1.0) = sin(Ο€ x)

The analytical solution is:

Ο•(x, y) = y sin(Ο€ x)

Reference: Olver PJ (2020) - page 210 - Introduction to Partial Differential Equations, Springer

See the code test_pde_poisson_2.rs

The results are illustrated below:

Poisson equation 2

Poisson equation 3

Approximate (with the Finite Differences Method, FDM) the solution of

βˆ‚Β²Ο•     βˆ‚Β²Ο•
β€”β€”β€”  +  β€”β€”β€” =  source(x, y)
βˆ‚xΒ²     βˆ‚yΒ²

on a (1.0 Γ— 1.0) square with homogeneous essential boundary conditions

The source term is given by (for a manufactured solution):

source(x, y) = 14yΒ³ - (16 - 12x) yΒ² - (-42xΒ² + 54x - 2) y + 4xΒ³ - 16xΒ² + 12x

The analytical solution is:

Ο•(x, y) = x (1 - x) y (1 - y) (1 + 2x + 7y)

See the code test_pde_poisson_3.rs

The results are illustrated below:

Poisson equation 3

Dependencies

~2.7–4.5MB
~75K SLoC