Time integrators

The ponio library provides different category of numerical integrators:

• Runge-Kutta methods based on Butcher tableaus,

• Extended stability methods,

• Splitting methods,

• IMEX methods with extended stability method.

For all examples in this section we will solve the Curtiss-Hirschfelder equation

\[\dot{y} = k(\cos(t) - y)\]

with stiff parameter \(k=50\), initial condition \(y(0) = y_0 = 2\) and time \(t\in[0, 4]\) with an initial time step \(\Delta t=0.05\).

    double const k   = 50.;
    double const y_0 = 2.;

    ponio::time_span<double> t_span = { 0., 4.0 };
    double dt                       = 0.05;

We will solve this problem with different kind of methods just to present how to use them.

Runge-Kutta methods from Butcher tableau

A Runge-Kutta method is defined in ponio by its Butcher tableau

\[\begin{split}\begin{array}{c|c} c & A \\ \hline & b^\top \end{array}\end{split}\]

which is stored as a JSON file in database folder.

Runge-Kutta methods are useful to solve a problem like

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

The method, with the previous Butcher tableau, reads as

\[\begin{split}\begin{aligned} u^{(i)} &= u^n + \Delta t \sum_j a_{ij} k_j, \quad i = 1, \dots, s \\ k_i &= f(t^n + c_i\Delta t, u^{(i)}) \\ u^{n+1} &= u^n + \Delta t \sum_i b_i k_i \end{aligned}\end{split}\]

Explicit methods

When the matrix \(A\) is strictly lower triangular, the Runge-Kutta method is called explicit. The only think you need to provide to solve a problem with this kind of method is the function \(f\).

See also

See the list of explicit Runge-Kutta methods in ponio.

For an explicit method, we can only define a the function \(f\) as a lambda as

auto f = [=]( double t, double y, double& dy )
{
    dy = k * ( std::cos( t ) - y );
};

And next call the ponio::solve() function with

ponio::solve( f, ponio::runge_kutta::rk_33(), y_0, t_span, dt, "ch_erk.txt"_fobs );

Embedded methods

Some Butcher tableaus provide an extra \(b\) line:

\[\begin{split}\begin{array}{c|c} c & A \\ \hline & b^\top \\ \hline & \hat{b}^\top \end{array}\end{split}\]

The Butcher tableau reads as

\[\begin{split}\begin{aligned} u^{(i)} &= u^n + \Delta t \sum_j a_{ij} k_j, \quad i = 1, \dots, s \\ k_i &= f(t^n + c_i\Delta t, u^{(i)}) \\ u^{n+1} &= u^n + \Delta t \sum_i b_i k_i \\ \tilde{u}^{n+1} &= u^n + \Delta t \sum_i \hat{b}_i k_i \end{aligned}\end{split}\]

with \(u^{n+1}\) and \(\tilde{u}^{n+1}\) two estimate of solution at time \(t^{n+1}\) of different order, we can compute a normalized error estimate with the following formula, with a given absolute tolerance \(a_{tol}\) respectively relative tolerance \(r_{tol}\)

\[e = \sqrt{ \frac{1}{N}\sum_i \left( \frac{|u^{n+1}_i - \tilde{u}^{n+1}_i|}{a_{tol} + r_{tol} \max( |u^n_i|, |u^{n+1}_i| )} \right)^2 }\]

This function is defined in ponio::detail::error_estimate(). Now compare the normalized error to \(1\) to accept or not the step, and compute a new time step

\[\begin{split}\begin{aligned} \Delta t_\text{new} &= 0.9 \sqrt[n]{ \frac{tol}{e} } \Delta t \\ \Delta t &= \min\left( \max\left( \Delta t_\text{new}, 0.2\Delta t \right) , 5\Delta t \right) \end{aligned}\end{split}\]

Most common embedded Runge-Kutta methods come from [DP80] and [PD81], that why tey are sometime call Dormand-Prince methods.

See also

See the list of embedded Runge-Kutta methods in ponio.

The most common embedded methods are explicit methods, so we can only define a the function \(f\) as a lambda as

auto f = [=]( double t, double y, double& dy )
{
    dy = k * ( std::cos( t ) - y );
};

And next call the ponio::solve() function with

ponio::solve( f, ponio::runge_kutta::rk54_6m().abs_tol( 1e-6 ).rel_tol( 1e-4 ), y_0, t_span, dt, "ch_dp.txt"_fobs );

Diagonal implicit methods

When the matrix \(A\) is lower triangular with a diagonal, the Runge-Kutta method is called diagonal implicit (or DIRK). You have to provide a Jacobian function that returns the Jacobian matrix in point \((t, u)\) (see ponio::implicit_problem). You can also provide an operator base definition (see ponio::implicit_operator_problem).

See also

See the list of diagonal implicit Runge-Kutta methods in ponio.

For an diagonal implicit Runge-Kutta method, we need the function \(f\) and also its Jacobian function \(\mathrm{d}f\), and store them into a ponio::implicit_problem.

auto f = [=]( double t, double y, double& dy )
{
    dy = k * ( std::cos( t ) - y );
};
auto df = [=]( double t, double /* y */ )
{
    return -k;
};

auto pb = ponio::make_implicit_problem( f, df );

And next call the ponio::solve() function with

ponio::solve( pb, ponio::runge_kutta::dirk34(), y_0, t_span, dt, "ch_dirk.txt"_fobs );

Lawson methods

Lawson methods was initially propose into [Law67]. Lawson methods are build to solve a problem with a linear and nonlinear part, to solve exactly the problem when the nonlinear part goes to zero. This class of problem can be write as

\[\dot{u} = L u + N(t, u)\]

First, we introduce the change of variable

\[v(t) = e^{-Lt}u(t)\]

which yields the equation

\[\dot{v}(t) = -Le^{-Lt}u(t) + e^{-Lt}\dot{u}(t)\]

which can be rewrite in therm of \(v\) as

\[\dot{v}(t) = e^{-Lt}N(t, e^{-Lt}v) = \tilde{N}(t, v)\]

We solve this equation with a classical Runge-Kutta method RK(\(s\), \(n\)) with \(s\) stages and of order \(n\). We rewrite the scheme in terms of the variable \(u\), which yields the following scheme

\[\begin{split}\begin{aligned} u^{(i)} &= u^n + \Delta t \sum_j a_{ij} k_j, \quad i = 1, \dots, s \\ k_i &= e^{-c_i\Delta t L} N(t^n + c_i\Delta t, e^{c_i\Delta t L} u^{(i)}) \\ u^{n+1} &= e^{\Delta t L}\left( u^n + \Delta t \sum_i b_i k_i \right) \end{aligned}\end{split}\]

In ponio, Lawson methods have the same name of the underlying explicit Runge-Kutta method prefixed by l.

See also

See the list of Lawson Runge-Kutta methods in ponio.

For a problem split into a linear and nonlinear part, we need to define linear part as a scalar or a matrix and the nonlinear part as a function, and store them into a ponio::lawson_problem.

double L = -k;
auto N   = [=]( double t, double y, double& dy )
{
    dy = k * std::cos( t );
};

auto pb = ponio::make_lawson_problem( L, N );

Because sometime we want to define exponential function with a Pade approximant, we need to specify an exponential function.

auto my_exp = []( double x )
{
    return std::exp( x );
};

And next call the ponio::solve() function with

ponio::solve( pb, ponio::runge_kutta::lrk_33( my_exp ), y_0, t_span, dt, "ch_lrk.txt"_fobs );

Exponential Runge-Kutta methods

Like Lawson methods, exponential Runge-Kutta methods are build to solve a problem with a linear and non-linear part, to solve exactly the problem when the nonlinear part goes to zero. This class of problem can be write as

\[\dot{u} = L u + N(t, u)\]

We solve this on a time step between \(0\) and \(\Delta t\)

\[u(t^n + \Delta t) = e^{\Delta t L} U + \int_0^{\Delta t} e^{(\Delta t - s)L}N(t^n + s, u(t^n+s)) \,\mathrm{d}s\]

Interpolation of the integral yields to build a custom Runge-Kutta method which reads as

\[\begin{split}\begin{aligned} u^{(i)} &= u^n + \Delta t\sum_j a_{ij}(\Delta t L)\cdot ( k_j + Lu^n ), \quad i = 1, \dots, s \\ k_i &= N(t^n + c_i\Delta t, u^{(i)}) \\ u^{n+1} &= u^n + \Delta t \sum_i b_i(\Delta t L)\cdot( k_i + Lu^n) \end{aligned}\end{split}\]

Note

Matrix \(A\) and vector \(b\) could contain some functions defined by

\[\varphi_\ell(z) = \frac{e^z - \sum_{k=0}^{\ell-1} \frac{1}{k!}z^k }{z^\ell}\]

and we use the notations \(\varphi_\ell = \varphi_\ell(\Delta t L)\) and \(\varphi_{\ell,j} = \varphi_\ell(c_j \Delta t L)\).

See also

See the list of exponential Runge-Kutta methods in ponio.

For a problem split into a linear and nonlinear part, we need to define linear part as a scalar or a matrix and the nonlinear part as a function, and store them into a ponio::lawson_problem.

double L = -k;
auto N   = [=]( double t, double y, double& dy )
{
    dy = k * std::cos( t );
};

auto pb = ponio::make_lawson_problem( L, N );

And next call the ponio::solve() function with

ponio::solve( pb, ponio::runge_kutta::exprk22(), y_0, t_span, dt, "ch_exprk.txt"_fobs );

Additive Runge-Kutta methods

An additive Runge-Kutta method is defined in ponio by its pair of Butcher tableaus

\[\begin{split}\begin{array}{c|c} c & A \\ \hline & b^\top \end{array} \quad \begin{array}{c|c} \tilde{c} & \tilde{A} \\ \hline & \tilde{b}^\top \end{array}\end{split}\]

which is stored as a JSON file in database folder, and where \((A, b, c)\) is an explicit Runge-Kutta method and \((\tilde{A}, \tilde{b}, \tilde{c})\) is a diagonal implicit Runge-Kutta method.

Additive Runge-Kutta methods are useful to solve a problem like

\[\dot{u} = f_e(t, u) + f_i(t,u)\]

The method, with the previous pair of Butcher tableaus, reads as

\[\begin{split}\begin{aligned} u^{(i)} &= u^n + \Delta t \sum_j a_{ij} k_j + \Delta t \sum_\ell \tilde{a}_{i\ell} \tilde{k}_\ell, \quad i = 1, \dots, s \\ k_i &= f_e(t^n + c_i\Delta t, u^{(i)}) \\ \tilde{k}_i &= f_i(t^n + \tilde{c}_i\Delta t, u^{(i)}) \\ u^{n+1} &= u^n + \Delta t \sum_i b_i k_i + \Delta t \sum_j \tilde{b}_j k_j \end{aligned}\end{split}\]

See also

See the list of additive Runge-Kutta methods in ponio.

For a problem split into a part which will be solved explicitly and a part which will be solved implicitly (you should provide the Jacobian matrix), and store them into a ponio::imex_problem:

auto f_e = [=]( double t, double y, double& dy )
{
    dy = k * std::cos( t );
};
auto f_i = [=]( double t, double y, double& dy )
{
    dy = -k * y;
};
auto df_i = [=]( double t, double y )
{
    return -k;
};

auto pb = ponio::make_imex_jacobian_problem( f_e, f_i, df_i );

And next call the ponio::solve() function with

ponio::solve( pb, ponio::runge_kutta::imex_rk36_spi2(), y_0, t_span, dt, "ch_ark.txt"_fobs );

---

Extended stability methods

Some problems, like heat equation, require methods stabilized on the negative real axis. The ponio library provides a Runge-Kutta Chebyshev method of order 2, ROCK2 method (of order 2), ROCK4 method (of order 4) and a Runge-Kutta Legendre method of order 1 or 2.

See also

See the list of extended stability Runge-Kutta methods in ponio.

Runge-Kutta Chebyshev method

The algorithm of RKC2 is the following:

\[\begin{split}\begin{aligned} y_0 &= f(t^n, y^n) \\ y_1 &= y_0 + \tilde{\mu}_1\Delta t f(t^n, y^n) \\ y_j &= (1 - \mu_j - \nu_j)y^n + \mu_j y_{j-1} + \nu_j y_{j-2} \\ &~~~~~~ + \tilde{\mu}_j \Delta t f(t^n + c_j\Delta t, y_{j-1}) + \tilde{\gamma}_j\Delta t f(t^n, y^n), \quad j=2,\dots, s \end{aligned}\end{split}\]

with coefficients given by

\[\tilde{\mu}_1 = b_1 \omega_1\]

\[\mu_j = \frac{2b_j}{b_{j-1}}\omega_0, \quad \nu_j = -\frac{b_j}{b_{j-2}}, \quad \tilde{\mu}_j = \frac{2b_j}{b_{j-1}}\omega_1\]

\[\tilde{\gamma}_j = -(1 - b_{j-1}T_{j-1}(\omega_0))\]

where

\[b_0 = b_2, \quad b_1 = \frac{1}{\omega_0}, \quad b_j = \frac{T_j''(\omega_0)}{(T_j'(\omega_0))^2},\ j=2,\dots, s\]

\[c_0 = 0, \quad c_1 = c_2, \quad c_j = \frac{T_s'(\omega_0)}{T_s''(\omega_0)}\frac{T_j''(\omega_0)}{T_j'(\omega_0)}, \quad c_s = 1\]

and

\[\omega_0 = 1 + \frac{\epsilon}{s^2},\quad \omega_1 = \frac{T_s'(\omega_0)}{T_s''(\omega_0)},\quad \epsilon \approx \frac{2}{13}\]

and where \(T_j(x)\) is the Chebyshev polynomial.

Runge-Kutta Chebyshev method is an explicit method, so we only need to define \(f\) function as

auto f = [=]( double t, double y, double& dy )
{
    dy = k * ( std::cos( t ) - y );
};

We need to specify how many stages we want (we choose 5), and next call the ponio::solve() function with

ponio::solve( f, ponio::runge_kutta::explicit_rkc2<5>(), y_0, t_span, dt, "ch_rkc.txt"_fobs );

ROCK2 method

We write the method ROCK2 presented in [AM01]. The algorithm of ROCK2 method is the following:

\[\begin{split}\begin{aligned} y_0 &= y^n \\ y_1 &= y^n + \Delta t \mu_1 f(y^n) \\ y_j &= \Delta t \mu_j f(y_{j-1}) - \nu_j y_{j-1} - \kappa_j y_{j-2}, \quad j=2,\dots, s-2\\ y_{s-1} &= y_{s-2} + \Delta t \sigma f(u_{s-2}) \\ y_{s}^\star &= y_{s-1} + \Delta t \sigma f(y_{s-1}) \\ y^{n+1} &= y_s^\star - \Delta t \sigma (1 - \frac{\tau}{\sigma})( f(y_{s-1}) - f(y_{s-2}) ) \end{aligned}\end{split}\]

where \(\mu_j\), \(\nu_j\) and \(\kappa_j\) coefficients coming from a minimization problem.

ROCK2 method is an explicit method, so we only need to define \(f\) function as

auto f = [=]( double t, double y, double& dy )
{
    dy = k * ( std::cos( t ) - y );
};

Next call the ponio::solve() function with

ponio::solve( f, ponio::runge_kutta::rock::rock2(), y_0, t_span, dt, "ch_rock2.txt"_fobs );

Hint

The ROCK2 method computes dynamically the number of stages with power method to estimate the spectral radius of your operator, but you can provide an other estimator of spectral radius if needed.

ROCK4 method

We write the method ROCK2 presented in [Abd02]. The algorithm of ROCK4 method is the following:

\[\begin{split}\begin{aligned} y_0 &= y^n \\ y_1 &= y^n + \Delta t \mu_1 f(y^n) \\ y_j &= \Delta t \mu_j f(y_{j-1}) - \nu_j y_{j-1} - \kappa_j y_{j-2}, \quad j=2,\dots, s-4 \\ y_{s-3} &= y_{s-4} + a_{21} \Delta t f(u_{s-4}) \\ y_{s-2} &= y_{s-4} + a_{31} \Delta t f(u_{s-4}) + a_{32} \Delta t f(y_{s-3}) \\ y_{s-1} &= y_{s-4} + a_{41} \Delta t f(u_{s-4}) + a_{42} \Delta t f(y_{s-3}) + a_{43} \Delta t f(y_{s-2}) \\ y^{n+1} &= y_{s-4} + b_1 \Delta t f(y_{s-4}) + b_2 \Delta t f(y_{s-3}) + b_3 \Delta t f(y_{s-2}) + b_4 \Delta t f(y_{s-1}) \end{aligned}\end{split}\]

where \(\mu_j\), \(\nu_j\) and \(\kappa_j\) coefficients coming from a minimization problem, and \(a_{ij}\), \(b_i\) coming from an order 4 method build with \(y_{s-4}\) as initial condition.

ROCK4 method is an explicit method, so we only need to define \(f\) function as

auto f = [=]( double t, double y, double& dy )
{
    dy = k * ( std::cos( t ) - y );
};

Next call the ponio::solve() function with

ponio::solve( f, ponio::runge_kutta::rock::rock2(), y_0, t_span, dt, "ch_rock4.txt"_fobs );

Hint

The ROCK4 method computes dynamically the number of stages with power method to estimate the spectral radius of your operator, but you can provide an other estimator of spectral radius if needed.

Runge-Kutta Legendre method

An other way to get a stabilized Runge-Kutta method is to use Legendre polynomials, we follow presentation in [MBA14].

The algorithm of RKL1 is the following:

\[\begin{split}\begin{aligned} y^{(0)} &= y^n \\ y^{(1)} &= y^n + \tilde{\mu}_1\Delta t f(t^n, y^n) \\ y^{(j)} &= \mu_jy^{(j-1)} + \nu_j y^{(j-2)} + \tilde{\mu}_j\Delta t f(t^n, y^{(j-1)}), \quad j=2,\dots, s \\ y^{(n+1)} &= y^{(s)} \end{aligned}\end{split}\]

with coefficients given by

\[\mu_j = \frac{2j-1}{j}, \qquad \nu_j = \frac{1-j}{j}\]

\[\tilde{\mu}_j = \frac{2j-1}{j}\frac{2}{s^2 + s}\]

where \(s\) is the number of stages of the method.

Runge-Kutta Legendre first-order method is an explicit method, so we only need to define \(f\) function as

auto f = [=]( double t, double y, double& dy )
{
    dy = k * ( std::cos( t ) - y );
};

We need to specify how many stages we want (we choose 5), and next call the ponio::solve() function with

ponio::solve( f, ponio::runge_kutta::explicit_rkl1<5>(), y_0, t_span, dt, "ch_rkl1.txt"_fobs );

The algorithm of RKL2 is the following:

\[\begin{split}\begin{aligned} y^{(0)} &= y^n \\ y^{(1)} &= y^{(0)} + \tilde{\mu}_1\Delta t f(t^n, y^{(0)}) \\ y^{(j)} &= \mu_jy^{(j-1)} + \nu_j y^{(j-2)} + (1-\mu_j-\nu_j)y^{(0)} + \tilde{\mu}_j\Delta t f(t^n, y^{(j-1)}) + \tilde{\gamma}_j\Delta tf(t^n, y^{(0)}), \quad j=2,\dots, s \\ y^{(n+1)} &= y^{(s)} \end{aligned}\end{split}\]

with coefficients given by

\[\mu_j = \frac{2j-1}{j}\frac{b_j}{b_{j-1}}, \qquad \nu_j = -\frac{j-1}{j}\frac{b_j}{b_{j-2}}\]

\[\tilde{\mu}_j = \mu_j w_1,\ 1<j \qquad \tilde{\mu}_1 = b_1 w_1\]

\[\tilde{\gamma}_j = -a_{j-1}\tilde{\mu}_j\]

where

\[b_0 = b_1 = b_2 = \frac{1}{3},\qquad b_j = \frac{j^2 + j - 2}{2j(j+1)},\ 2\leq j\]

and

\[a_j = 1 - b_j, \qquad w_1 = \frac{4}{s^2 + s - 2}\]

where \(s\) is the number of stages of the method

Runge-Kutta Legendre second-order method is an explicit method, so we only need to define \(f\) function as

auto f = [=]( double t, double y, double& dy )
{
    dy = k * ( std::cos( t ) - y );
};

We need to specify how many stages we want (we choose 5), and next call the ponio::solve() function with

ponio::solve( f, ponio::runge_kutta::explicit_rkl2<5>(), y_0, t_span, dt, "ch_rkl2.txt"_fobs );

---

Splitting methods

Some methods to solve a problem of the form:

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

with the initial condition \(u(t=0)=u_0\). We note with \(\phi_{\tau}^{[f_i]}(t^n,\tilde{u}^n)\) the solution at time \(t^n+\tau\) of the subproblem \(\dot{u}=f_i(t,u)\) with the initial condition \(u(t^n)=\tilde{u}^n\).

See also

See the list of splitting methods in ponio.

Lie Splitting method

In Lie splitting method, the solution is computed as:

\[u^{n+1} = \mathcal{L}^{\Delta t}(t^n, u^n) = \phi_{\Delta t}^{[f_1]}\circ \cdots \circ \phi_{\Delta t}^{[f_k]} (t^n,u^n)\]

For splitting methods in ponio, we first need to define all sub-problem (here only \(f_1\) and \(f_2\)) and store them into a ponio::problem

auto f1 = [=]( double t, double y, double& dy )
{
    dy = k * std::cos( t );
};
auto f2 = [=]( double t, double y, double& dy )
{
    dy = -k * y;
};

auto pb = ponio::make_problem( f1, f2 );

We also define the Lie splitting method and the method to solve each substep. In this example we use two RK(3,3) methods, and we specify the time step for subcycle

auto lie = ponio::splitting::make_lie_tuple( std::make_pair( ponio::runge_kutta::rk_33(), 0.5 * dt ),
    std::make_pair( ponio::runge_kutta::rk_33(), 0.5 * dt ) );

Next call the ponio::solve() function with

ponio::solve( pb, lie, y_0, t_span, dt, "ch_lie.txt"_fobs );

Strang Splitting method

In Strang splitting method, the solution is computed as:

\[u^{n+1} = \mathcal{S}^{\Delta t}(t^n, u^n) = \phi_{\frac{\Delta t}{2}}^{[f_1]}\circ \cdots \circ \phi_{\frac{\Delta t}{2}}^{[f_{k-1}]} \circ \phi_{\Delta t}^{[f_k]} \circ \phi_{\frac{\Delta t}{2}}^{[f_{k-1}]}\circ\cdots\circ \phi_{\frac{\Delta t}{2}}^{[f_1]} (t^n,u^n)\]

For splitting methods in ponio, we first need to define all sub-problem (here only \(f_1\) and \(f_2\)) and store them into a ponio::problem

auto f1 = [=]( double t, double y, double& dy )
{
    dy = k * std::cos( t );
};
auto f2 = [=]( double t, double y, double& dy )
{
    dy = -k * y;
};

auto pb = ponio::make_problem( f1, f2 );

We also define the Strang splitting method and the method to solve each substep. In this example we use two RK(3,3) methods, and we specify the time step for subcycle

auto strang = ponio::splitting::make_strang_tuple( std::make_pair( ponio::runge_kutta::rk_33(), 0.5 * dt ),
    std::make_pair( ponio::runge_kutta::rk_33(), 0.5 * dt ) );

Next call the ponio::solve() function with

ponio::solve( pb, strang, y_0, t_span, dt, "ch_strang.txt"_fobs );

Adaptive time step Strang splitting method

In Strang splitting method, the solution is computed as:

\[u^{n+1} = \mathcal{S}^{\Delta t}(t^n, u^n) = \phi_{\frac{\Delta t}{2}}^{[f_1]}\circ \cdots \circ \phi_{\frac{\Delta t}{2}}^{[f_{k-1}]} \circ \phi_{\Delta t}^{[f_k]} \circ \phi_{\frac{\Delta t}{2}}^{[f_{k-1}]}\circ\cdots\circ \phi_{\frac{\Delta t}{2}}^{[f_1]} (t^n,u^n)\]

The approximation of order 1 is computed with a shifted Strang splitting method:

\[\tilde{u}^{n+1} = \mathcal{S}_{\delta}^{\Delta t}(t^n, u^n) = \phi_{(\frac{1}{2}-\delta)\Delta t}^{[f_1]}\circ\phi_{\frac{\Delta t}{2}}^{[f_2]}\circ \cdots \circ \phi_{\frac{\Delta t}{2}}^{[f_{k-1}]} \circ \phi_{\Delta t}^{[f_k]} \circ \phi_{\frac{\Delta t}{2}}^{[f_{k-1}]}\circ\cdots\circ\phi_{\frac{\Delta t}{2}}^{[f_2]}\circ \phi_{(\frac{1}{2}+\delta)\Delta t}^{[f_1]} (t^n,u^n)\]

with \(\delta\in[-1/2, 0)\cup(0,1/2]\)

Note

Only the first (and last) step is shifted by \(\delta\Delta t\).

The difficulty in adaptive time step Strang splitting method is to find a good shift coefficient \(\delta\) for a given tolerance \(\eta\). The complete strategy of an iteration provides in [DDM15] is

1. With a given tolerance \(\eta\), time step \(\Delta t\) and shift \(\delta\) compute estimates

   \[\begin{split}\begin{aligned} u^{n+1} &= \mathcal{S}^{\Delta t}(t^n, u^n) \\ \tilde{u}^{n+1} &= \mathcal{S}_{\delta}^{\Delta t}(t^n, u^n) \end{aligned}\end{split}\]

2. Compute new time step from local error

   \[\Delta t^{new} = s \Delta t \sqrt{\frac{\eta}{\| u^{n+1} - \tilde{u}^{n+1} \|}}\]

3. If local error is lower than tolerance: \(\| u^{n+1} - \tilde{u}^{n+1} \| < \eta\), iteration is accepted, else compute a new iteration with the new time step.

4. Every \(N_\delta\) iterations, or if \(\Delta t\not\in[\beta\Delta t^\star, \gamma\Delta t^\star]\) (or if first iteration and \(\Delta t^\star\) is not computed), compute new \(\Delta t^\star\) from current \(\delta\) (given by info() returned value which has a data member variable) and \(C_0\) (given by ponio::splitting::strang::adaptive_strang::lipschitz_constant_estimate() method) with:

   \[\Delta t^\star \approx \frac{\delta C_\delta}{C_0},\qquad \delta C_\delta \Delta t^2 = \| \mathcal{S}^{\Delta t}(t^n, u^n) - \mathcal{S}_\delta^{\Delta t}(t^n, u^n) \|\]

5. If \(\Delta t \not\in [\beta\Delta t^\star, \gamma\Delta t^\star]\), assume \(\Delta t^\star=\Delta t\) and compute a new \(\delta\) with:

   \[\Delta t^\star \approx \frac{\delta C_\delta}{C_0}\]

---

IMEX methods with extended stability method

The PIROCK method is introduce in [AV13], the complete scheme is a IMEX scheme that allows for solving an equation of the following form (with 3 operators):

\[\dot{u} = F_R(u) + F_D(u) + F_A(u)\]

with \(F_R\) a reaction operator (solved implicitly), \(F_D\) a diffusion operator (solved explicitly by a modified ROCK2 method) and \(F_A\) an advection operator (solved with an explicit RK3 method). The first implementation of PIROCK scheme in ponio makes to solve only reaction-diffusion problem (e.g. \(F_A = 0\)), the second one to solve a complet reaction-diffusion-advection problem.

The method is divided into 5 steps:

1. Diffusion stabilization procedure (modified ROCK2 method)

   \[\begin{split}\begin{aligned} u^{(0)} &= u^n \\ u^{(1)} &= u^n + \alpha \mu_1 \Delta t F_D(u^n) \\ u^{(j)} &= \alpha \mu_j \Delta t F_D(u^{(j-1)}) - \nu_j u^{(j-1)} - \kappa_j u^{(j-2)},\quad j=2,\dots, s-2+\ell \end{aligned}\end{split}\]

2. Finishing procedure for diffusion

   \[\begin{split}\begin{aligned} u^{*(s-1)} &= u^{(s-2)} + \sigma_\alpha \Delta t F_D(u^{(s-2)}) \\ u^{*(s)} &= u^{*(s-1)} + \sigma_\alpha \Delta t F_D(u^{*(s-1)}) \end{aligned}\end{split}\]

3. Starting value for advection-reaction

   \[U = u^{(s-2+\ell)}\]

4. Finishing procedure for advection-reaction coupling

   \[\begin{split}\begin{aligned} u^{(s+1)} &= u^{(s-2+\ell)} + \gamma \Delta t F_R(u^{(s+1)}) \\ u^{(s+2)} &= u^{(s-2+\ell)} + \beta \Delta t F_D(u^{(s+1)}) + \Delta t F_A(u^{(s+1)}) + (1-2\gamma)\Delta t F_R(u^{(s+1)}) + \gamma \Delta t F_R(u^{(s+2)}) \\ u^{(s+3)} &= u^{(s-2+\ell)} + (1-2\gamma)\Delta t F_A(u^{(s+1)}) + (1-\gamma)\Delta t F_R(u^{(s+1)}) \\ u^{(s+4)} &= u^{(s-2+\ell)} + \frac{1}{3}\Delta t F_A(u^{(s+1)}) \\ u^{(s+5)} &= u^{(s-2+\ell)} + \frac{2\beta}{3} \Delta t F_D(u^{(s+1)}) + \frac{2}{3}\Delta t J_R^{-1} F_A(u^{(s+4)}) + \left(\frac{2}{3} - \gamma\right) \Delta t F_R(u^{(s+1)}) + \frac{2\gamma}{3}\Delta t F_R(u^{(s+2)}) \end{aligned}\end{split}\]

5. Computation of the integration step

   \[\begin{split}\begin{aligned} u^{n+1} =& u^{*(s)} \\ & - \sigma_\alpha\left( 1 - \frac{\tau_\alpha}{\sigma_\alpha^2}\right)\Delta t (F_D(u^{*(s-1)}) - F_D(u^{(s-2)})) \\ & + \frac{1}{4}\Delta t F_A(u^{(s+1)}) + \frac{3}{4}\Delta t F_A(u^{(s+5)}) \\ & + \frac{1}{2}\Delta t F_R(u^{(s+1)}) + \frac{1}{2}\Delta t F_R(u^{(s+2)}) \\ & + \frac{J_R^{-\ell}}{2-4\gamma}\Delta t (F_D(u^{(s+3)}) - F_D(u^{(s+1)})) \end{aligned}\end{split}\]

   where \(\mu_j\), \(\nu_j\), \(\kappa_j\) are the same coefficients as for the standard ROCK2 method and:

   • \(\gamma = 1-\frac{\sqrt{2}}{2}\)

   • \(\beta = 1-2\alpha P'_{s-2+\ell}(0)\) where \(P'_{s-2+\ell}\) is the stability polynomial of the underlying ROCK2 method

   • \(J_R = I - \gamma \Delta t \frac{\partial F_R}{\partial u}(u^{(s-2+\ell)})\)

   • \(\sigma_\alpha = \frac{1-\alpha}{2} + \alpha \sigma\), where \(\sigma\) coming from the underlying ROCK2 method

   • \(\tau_\alpha = \frac{(\alpha-1)^2}{2} + 2\alpha(1-\alpha)\sigma + \alpha^2\tau\), where \(\sigma\) and \(\tau\) coming from the underlying ROCK2 method

For adaptive time step method, need to compute one error by operator and take the maximum value

\[\begin{split}\begin{aligned} err_D &= \sigma_\alpha \Delta t \left( 1 - \frac{\tau_\alpha}{\sigma_\alpha^2} \right) ( F_D( u^{*(s-1)} - u^{(s-2)} )) \\ err_R &= \Delta t J_R^{-1}\left( \frac{1}{6}F_R(u^{(s+1)}) - \frac{1}{6}F_R(u^{(s+2)}) \right) \\ err_A &= -\frac{3}{20}\Delta t F_A(u^{(s+1)}) + \frac{3}{10}\Delta t F_A(u^{(s+4)}) - \frac{3}{20}\Delta t F_A(u^{(s+5)}) \end{aligned}\end{split}\]

The error is computed with the maximum:

\[err = \max\left( \|err_D\|, \|err_R\|, \|err_A\|^{2/3} \right)\]

with the following error norm:

\[\|\cdot\| : err \mapsto \sum_i \left( \frac{err_i}{a_{tol} + \max(y^{n+1}_i, y^n_i) } \right)^2\]

where \(y^n\) and \(y^{n+1}\) are estimation of the solution respectively at time \(t^n\) and after an iteration \(t^{n+1}=t^n+\Delta t\), and \(a_{tol}\) and \(r_{tol}\) respectively absolute and relative tolerances. This norm give a normalized error, so its compare to \(1\) and new time step is computed as:

\[\Delta t_{new} = \sqrt{\frac{1}{err}} \Delta t\]

Note

As in a lot of adaptive time step method, the new time step stays in \(\Delta t_{new} \in [0.5\Delta t, 2 \Delta t]\) set and there is a safety factor, so the new time step is multiply by \(0.8\).

Still keep two free parameters \(\ell\) and \(\alpha\) given free to user, ponio provides two choices for this parameters:

• \(\ell=2\) and \(\alpha=1\), in this case, if \(F_A=0\) and \(F_R=0\) we have the standard ROCK2 method;

• \(\ell=1\) and \(\alpha = \frac{1}{2P'_{s-2+\ell}(0)}\), so \(\beta=0\) to minimized computation cost.

See also

See the list of IMEX methods with extended stability method in ponio.

For a PIROCK method to solve only reaction-diffusion problem, we need two functions \(F_D\) and \(F_R\) and the Jacobian function \(\partial_u F_R\), and store them into a ponio::imex_problem.

auto f1 = [=]( double t, double y, double& dy )
{
    dy = k * std::cos( t );
};
auto f2 = [=]( double t, double y, double& dy )
{
    dy = -k * y;
};
auto df2 = [=]( double t, double y )
{
    return -k;
};

auto pb = ponio::make_imex_jacobian_problem( f1, f2, df2 );

And next call the ponio::solve() function with

ponio::solve( pb, ponio::runge_kutta::pirock::pirock_a1(), y_0, t_span, dt, "ch_pirock.txt"_fobs );

---

Bibliography

[Abd02]

Assyr Abdulle. Fourth order chebyshev methods with recurrence relation. SIAM Journal on Scientific Computing, 23(6):2041–2054, 2002. doi:10.1137/S1064827500379549.

[AM01]

Assyr Abdulle and Alexei A. Medovikov. Second order chebyshev methods based on orthogonal polynomials. Numerische Mathematik, 2001. doi:10.1007/s002110100292.

[AV13]

Assyr Abdulle and Gilles Vilmart. Pirock: a swiss-knife partitioned implicit–explicit orthogonal runge–kutta chebyshev integrator for stiff diffusion–advection–reaction problems with or without noise. Journal of Computational Physics, 242:869–888, 2013. doi:10.1016/j.jcp.2013.02.009.

[DDM15]

Stéphane Descombes, Max Duarte, and Marc Massot. Operator splitting methods with error estimator and adaptive time-stepping. Application to the simulation of combustion phenomena. In Roland Glowinski, Stanley Osher, and Wotao Yin, editors, Splitting Methods in Communication, Imaging, Science, and Engineering, pages 1–13. Springer International Publishing, 2015. doi:10.1007/978-3-319-41589-5.

[DP80]

J.R. Dormand and P.J. Prince. A family of embedded runge-kutta formulae. Journal of Computational and Applied Mathematics, 6(1):19–26, 1980. doi:10.1016/0771-050X(80)90013-3.

[Law67]

J. Douglas Lawson. Generalized runge-kutta processes for stable systems with large lipschitz constants. SIAM Journal on Numerical Analysis, 4(3):372–380, 1967. doi:10.1137/0704033.

[MBA14]

Chad D. Meyer, Dinshaw S. Balsara, and Tariq D. Aslam. A stabilized runge–kutta–legendre method for explicit super-time-stepping of parabolic and mixed equations. Journal of Computational Physics, 257:594–626, 2014. doi:10.1016/j.jcp.2013.08.021.

[PD81]

P.J. Prince and J.R. Dormand. High order embedded runge-kutta formulae. Journal of Computational and Applied Mathematics, 7(1):67–75, 1981. doi:10.1016/0771-050X(81)90010-3.