ponio#

The library ponio is a collection of time integrators for solving differential equations written in C++

Link

Github repository

Language

C++

License

BSD-3-Clause

Documentation

documentation

Analysis

analysis of methods

Lorenz equations#

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

\[\begin{split} \begin{cases} \dot{x} &= \sigma (y - x) \\ \dot{y} &= \rho x - y - x z \\ \dot{z} &= x y - \beta z \end{cases} \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 ponio::runge_kutta::rk_44 method.

To solve a problem with ponio we first write our equation as a ODE of the form:

\[ \dot{u} = f(t, u) \]

and the user provide the function \(f\) as (a lambda function for the example):

  auto f = [](double t, auto const& u, state_t& du);

where u the current state of the function, t the current time and output is the result of \(f(t,u)\).

19    auto lorenz = [&]( double t, auto&& y, state_t& dy )
20    {
21        dy[0] = sigma * ( y[1] - y[0] );
22        dy[1] = y[0] * ( rho - y[2] ) - y[1];
23        dy[2] = y[0] * y[1] - beta * y[2];
24    };

After defined a method with ponio::runge_kutta::rk_44(), we can solve the problem between initial time and final time and give an observer which be call after each succeed time iteration

32    ponio::solve( lorenz, ponio::runge_kutta::rk_44(), y0, t_span, dt, ponio::observer::file_observer( "lorenz.txt" ) );

the ponio::observer::file_observer is an observer that save all states in a file.

For the complet example, see lorenz.cpp source file.

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 the initial condition given by a hat function:

\[\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 and a forward Euler method for the time discretization given by ponio::runge_kutta::euler method.

We define the up-wind scheme as:

50    auto upwind = [a, n_x, dx]( double t, auto&& y, state_t& dy )
51    {
52        dy[0] = -( std::max( a, 0. ) * ( y[0] - y[n_x - 1] ) + std::min( a, 0. ) * ( y[1] - y[0] ) ) / dx;
53
54        for ( std::size_t i = 1; i < n_x - 1; ++i )
55        {
56            dy[i] = -( std::max( a, 0. ) * ( y[i] - y[i - 1] ) + std::min( a, 0. ) * ( y[i + 1] - y[i] ) ) / dx;
57        }
58
59        dy[n_x - 1] = -( std::max( a, 0. ) * ( y[n_x - 1] - y[n_x - 2] ) + std::min( a, 0. ) * ( y[0] - y[n_x - 1] ) ) / dx;
60    };

The time loop is the same as for Lorenz equation.

For the complet example, see transport.cpp source file.

Arenstorf orbit#

The Arenstorf orbit problem is a classical problem to test adaptive time step methods:

\[\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}\]

We define this system as:

17    auto arenstorf = [=]( double t, auto&& y, state_t& dy )
18    {
19        double const y1 = y[0];
20        double const y2 = y[1];
21        double const y3 = y[2];
22        double const y4 = y[3];
23
24        double const r1 = sqrt( ( y1 + mu ) * ( y1 + mu ) + y2 * y2 );
25        double const r2 = sqrt( ( y1 - 1 + mu ) * ( y1 - 1 + mu ) + y2 * y2 );
26
27        dy[0] = y3;
28        dy[1] = y4;
29        dy[2] = y1 + 2 * y4 - ( 1 - mu ) * ( y1 + mu ) / ( r1 * r1 * r1 ) - mu * ( y1 - 1 + mu ) / ( r2 * r2 * r2 );
30        dy[3] = y2 - 2 * y3 - ( 1 - mu ) * y2 / ( r1 * r1 * r1 ) - mu * y2 / ( r2 * r2 * r2 );
31    };

We solve this example with given method ponio::runge_kutta::rk54_7m which is the method RK5(4) 7M in [DP80] (mainly call DOPRI5) and ponio::runge_kutta::rk87_13m which is the method RK8(7) 13M in [PD81] (mainly call DOPRI8).

The time loop is the same as for Lorenz equation, for rk54_7m method

39    ponio::solve( arenstorf, ponio::runge_kutta::rk54_7m( 1e-5 ), y0, t_span, dt, ponio::observer::file_observer( "arenstorf_rk54_7m.txt" ) );

and for rk87_13m method

40    ponio::solve( arenstorf, ponio::runge_kutta::rk87_13m( 1e-5 ), y0, t_span, dt, ponio::observer::file_observer( "arenstorf_rk87_13m.txt" ) );

For the complet example, see arenstorf.cpp source file.