### 13 stable releases

new 1.6.0 | May 21, 2024 |
---|---|

1.5.1 | May 16, 2024 |

1.4.0 | Apr 30, 2024 |

1.0.0 | Mar 30, 2024 |

0.9.2 | Mar 17, 2024 |

#**54** in Math

**1,339** downloads per month

**MIT**license

**3MB**

**40K**
SLoC

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

*This crate is part of Russell - Rust Scientific Library*

## Contents

- Introduction
- Installation
- π Examples

## 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:

for ODE systems and non-stiff problems using moderate tolerances`DoPri5`

for ODE systems and non-stiff problems using strict tolerances`DoPri8`

for ODE and DAE systems, possibly stiff, with moderate to strict tolerances`Radau5`

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

### Documentation

### References

- 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
- 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
- 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

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

`[``dependencies``]`
`russell_ode ``=` `"`*`"`

### Optional features

The following (Rust) features are available:

: Use Intel MKL instead of OpenBLAS`intel_mkl`

: Use a locally compiled version of SuiteSparse`local_suitesparse`

: Enable the MUMPS solver (locally compiled)`with_mumps`

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

. See also:`russell_ode`

### 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`.``init_mass_matrix``(`mass_nnz`,` symmetric`)``?``;`
system`.``mass_put``(``0``,` `0``,` `1.``0``)``?``;`
system`.``mass_put``(``0``,` `1``,` `1.``0``)``?``;`
system`.``mass_put``(``1``,` `0``,` `1.``0``)``?``;`
system`.``mass_put``(``1``,` `1``,` `-``1.``0``)``?``;`
system`.``mass_put``(``2``,` `2``,` `1.``0``)``?``;`
`//` 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:

#### Variable step sizes

This example solves the Brusselator ODE with variable step sizes for different tolerances. In this example,

. The global error is the difference between Russell's results and Mathematica's results obtained with high accuracy. The Mathematica code is:`tol = abs_tol = rel_tol`

`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:

#### 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 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

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.`second_book`

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

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

` β§ 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

.
However, here we assume periodic on `u`

and `u`

:`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

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

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

and `U`

as follows:`V`

`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

, the components of the Jacobian matrix can be "assembled" as follows:`Fβ := βYβ/βt`

` β§ β§ β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

and is implemented in brusselator_pde_radau5.rs. This code will generate a series of files, one for each (dense) output time with `npoint = 21`

`h_out ``=` `0.``5`

.The results can then be plotted with brusselator_pde_plot.rs

The results are shown below for the

field:`U`

And below for the

field:`V`

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

The computations with

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

The figure below shows the

(black dashed lines) and Mathematica (red solid lines) results for the `russell`

field:`U`

The figure below shows the

(black dashed lines) and Mathematica (red solid lines) results for the `russell`

field:`V`

#### Second book

For the problem in the second book, we run brusselator_pde_radau5_2nd.rs with

and `npoint = 129`

`h_out ``=` `1.``0`

.The results are shown below for the

field:`U`

And below for the

field:`V`

The code brusselator_pde_2nd_comparison.rs compares

results with Mathematica results.`russell`

The figure below shows the

(black dashed lines) and Mathematica (red solid lines) results for the `russell`

field:`U`

The figure below shows the

(black dashed lines) and Mathematica (red solid lines) results for the `russell`

field:`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:

### 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:

### 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:

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

### 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:

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:

### 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:

### PDE: discrete Laplacian operator in 2D

For convenience (e.g., in benchmarks),

implements a discrete Laplacian operator (2D) based on the Finite Differences Method.`russell_ode`

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:

#### 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 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 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: