Comparison#

This section is the comparison with some ODE solvers, we present the same examples writing with other ODE solvers to show choices of interface in ponio.

Lorenz equations#

We would like to solve the Lorenz equations, a classical chaotic system given by:

\[\begin{split}\begin{aligned} \dot{x} &= \sigma ( y - x ) \\ \dot{y} &= \rho x - y - x z \\ \dot{z} &= x y - \beta z \end{aligned}\end{split}\]

with parameter \(\sigma = 10\), \(\rho = 28\) and \(\beta = \frac{8}{3}\), and the initial state \((x_0, y_0, z_0) = (1, 1, 1)\). We solve it with a classical Runge-Kutta method of order 4, given by its Butcher tableau:

\[\begin{split}\begin{array}{c|cccc} 0 & & & & \\ \frac{1}{2} & \frac{1}{2} & & & \\ \frac{1}{2} & 0 & \frac{1}{2} & & \\ 1 & 0 & 0 & 1 & \\ \hline & \frac{1}{6} & \frac{1}{3} & \frac{1}{3} & \frac{1}{6} \end{array}\end{split}\]

This method is pretty much always in library with various name:

except in SciPy where we need to add this.

Transport equation#

In this example we would like to solve the following PDE:

\[\partial_t u + a \partial_x u = 0\]

with \(t>0\), on the torus \(x\in [0, 1)\), a velocity \(a=1\) and with the initial condition given by:

\[\begin{split}u(0, x) = \begin{cases} x - 0.25 & \text{if } x\in[0.25, 0.5[ \\ -x + 0.75 & \text{if } x\in[0.5, 0.75[ \\ 0 & \text{else} \end{cases}\end{split}\]

We choose a first order up-wind scheme to estimate the \(x\) derivative:

\[\begin{split}a \partial_x u(t^n, x_i) \approx \begin{cases} \frac{u^n_{i+1} - u^n_{i}}{\Delta x}, & \text{if } a \geq 0 \\ \frac{u^n_{i} - u^n_{i-1}}{\Delta x}, & \text{else} \end{cases}\end{split}\]

and the forward Euler method for time discretization:

\[u^{n+1} = u^n - \Delta t \textrm{D}_a(u^n)\]

where \(\textrm{D}_a\) is the approximation of partial derivative in \(x\) at velocity \(a\). The complet scheme, with \(a>0\), becomes:

\[u^{n+1}_i = u^n_i - a \frac{\Delta t}{\Delta t}( u^n_{i+1} - u^n_i )\]

In case of \(a=1\) and \(\Delta t = \Delta x\), the scheme gives the exact solution. The explicit Euler method (or forward Euler method) is present in:

  • Ascent with the name asc::Euler

  • odeint with the name euler

  • petsc with the name TSRK1FE

  • ponio with the name euler

In GSL we need to implement our own driver, and in SciPy to provide a Butcher tableau.

../_images/transport.png

Arenstorf orbit#

We would like to compare some adaptive time step methods with the Arenstorf orbit problem

\[\begin{split}\begin{cases} \ddot{x} &= x + 2\dot{y} - \frac{1-\mu}{r_1^3}(x+\mu) - \frac{\mu}{r_2^3}(x-1+\mu) \\ \ddot{y} &= y - 2\dot{x} - \frac{}{1-\mu}{r_1^3}y - \frac{\mu}{r_2^3}y \end{cases}\end{split}\]

with initial condition \((x,\dot{x},y,\dot{y})=(0.994, 0, 0, -2.001585106)\), \(r_1\) and \(r_2\) given by

\[r_1 = \sqrt{(x+\mu)^2 + y^2},\quad r_2 = \sqrt{(x-1+\mu)^2 + y^2}\]

and with parameter \(\mu = 0.012277471\).

First of all, we need to rewrite this problem into a first order derivative equation in time

\[\begin{split}\begin{pmatrix} y_1 \\ y_2 \\ y_3 \\ y_4 \end{pmatrix} = \begin{pmatrix} x \\ y \\ \dot{x} \\ \dot{y} \end{pmatrix}, \qquad \begin{cases} \dot{y}_1 = y_3 \\ \dot{y}_2 = y_4 \\ \dot{y}_3 = y_1 + 2y_4 - \frac{1-\mu}{r_1^3}(y_1 + \mu) - \frac{\mu}{r_2^3}(y_1-1+\mu) \\ \dot{y}_4 = y_2 - 2y_3 - \frac{1-\mu}{r_1^3}y_2 - \frac{\mu}{r_2^3}y_2 \\ \end{cases}\end{split}\]

This is a good example to use an adaptive time step method, because we need a small time step close to the Moon, elsewhere the constraint can be relaxed. For the sake of simplicity we use only native adaptive time step method in each library, it will be:

where DOPRI45, runge_kutta_dopri5, TSRK5DP, rk54_7m or RK45 are the same method, in [DP80] (RK5(4) 7M), see ponio::runge_kutta::rk54_7m_t for more information; and gsl_odeiv2_step_rk8pd, DOP853 and rk87_13m are the same method in [PD81] (RK8(7) 13M) ponio::runge_kutta::rk87_13m_t.

../_images/arenstorf.png