List of extended stabilized Runge-Kutta 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.

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.

template<std::size_t N_stages_, typename _value_t = double>
class explicit_rkc2#

define RKC2 with user defined number of stages

Template Parameters:

  • N_stages_ – number of stages

  • value_t – type of coefficients

Public Functions

inline explicit_rkc2(value_t eps = 2. / 13.)#

Construct a new explicit RKC2 algorithm.

Parameters:

eps – value of relaxation

template<typename problem_t, typename state_t, typename array_ki_t, std::size_t j>
inline void stage(Stage<j>, problem_t &f, value_t tn, state_t &yn, array_ki_t const &Yj, value_t dt, state_t &ui, state_t &yi)#

computes stage \(j\) of RKC2 method

\(y_j = (1 - \mu_j - \nu_j)y^n + \mu_j y_{j-1} + \nu_j y_{j-2} + \tilde{\mu}_j\Delta tf(t^n + c_{j-1}\Delta t, y_{j-1}) + \tilde{\gamma}_j\Delta t f(t^n, y^n)\)

Template Parameters:

  • problem_t – type of operator \(f\)

  • state_t – type of current state

  • array_ki_t – type of array of temporary stages

  • j – integer of stage

Parameters:

  • f – operator \(f\)

  • tn – current time

  • yn – current state

  • Yj – array of temporary stages

  • dt – current time step

  • ui – temporary step

  • yi – computed output step

template<typename problem_t, typename state_t, typename array_ki_t>
inline void stage(Stage<0>, problem_t &f, value_t tn, state_t &yn, array_ki_t const&, value_t, state_t&, state_t &yi)#

computes pseudo first stage of RKC2 method

\(y_0 = f(t^n, y^n)\)

Template Parameters:

  • problem_t – type of operator \(f\)

  • state_t – type of current state

  • array_ki_t – type of array of temporary stages

  • j – integer of stage

Parameters:

  • f – operator \(f\)

  • tn – current time

  • yn – current state

  • yi – computed output step

template<typename problem_t, typename state_t, typename array_ki_t>
inline void stage(Stage<1>, problem_t&, value_t, state_t &yn, array_ki_t const &Yj, value_t dt, state_t&, state_t &yi)#

computes first stage of RKC2 method

\(y_1 = y^n + \tilde{\mu}_1 \Delta t f(t^n, y^n)\)

Template Parameters:

  • problem_t – type of operator \(f\)

  • state_t – type of current state

  • array_ki_t – type of array of temporary stages

  • j – integer of stage

Parameters:

  • yn – current state

  • Yj – array of temporary stages

  • dt – current time step

  • yi – computed output step

template<typename problem_t, typename state_t, typename array_ki_t>
inline void stage(Stage<2>, problem_t &f, value_t tn, state_t &yn, array_ki_t const &Yj, value_t dt, state_t &ui, state_t &yi)#

computes second stage of RKC2 method

\(y_2 = (1-\mu_2 -\nu_2)y^n + \mu_2y_1 + \nu_2y^n + \tilde{\mu}_2\Delta t f(t^n + c_1\Delta t, y_1) + \tilde{\gamma}_2 \Delta t f(t^n, y^n)\)

Template Parameters:

  • problem_t – type of operator \(f\)

  • state_t – type of current state

  • array_ki_t – type of array of temporary stages

  • j – integer of stage

Parameters:

  • f – operator \(f\)

  • tn – current time

  • yn – current state

  • Yj – array of temporary stages

  • dt – current time step

  • ui – temporary step

  • yi – computed output step

inline auto &info()#

gets iteration_info object

inline auto const &info() const#

gets iteration_info object (constant version)

Public Static Functions

template<std::size_t J>
static inline value_t b(value_t const x)#

computes \(b_j = \) coefficients with following formula: \(b_0=b_2\), \(b_1=\frac{1}{\omega_0}\), \(b_j = \frac{T_j''(\omega_0)}{(T_j'(\omega_0))^2}\)

Template Parameters:

J – index \(j\)

Parameters:

x – value of \(\omega_0\)

Warning

In ponio library the number of stages of RKC2 is static, you can’t get dynamic number of stages for this method, look at ponio::runge_kutta::rock::rock2() or ponio::runge_kutta::rock::rock4() methods for dynamic number of stages (and optional adaptive time step method).

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.

template<typename eig_computer_t, bool _is_embedded = false, typename _value_t = double>
class rock2_impl#

define ROCK2 method

Template Parameters:

  • eig_computer_t – type of computer of maximal eigenvalue

  • _is_embedded – define if method is used as adaptive or constant time step method [default is false]

  • _value_t – type of coefficients

Public Functions

inline rock2_impl(eig_computer_t &&_eig_computer)#

Construct a new ROCK2 algorithm object.

Parameters:

_eig_computer – estimator of spectral radius

template<typename state_t>
inline auto error(state_t &&unp1, state_t &&un, state_t &&tmp)#

computes an error for adaptive time step method

Template Parameters:

state_t – type of current state

Parameters:

  • unp1 – estimation of solution at time \(t^{n+1} = t^n+\Delta t\)

  • un – solution at time \(t^n\)

  • tmp – \(-\Delta t \sigma(1 - \frac{\tau}{\sigma^2})(f(u_{s-1}) - f(u_{s-2}))\)

Returns:

auto estimation of error to compare to 1

template<typename problem_t, typename state_t, typename array_ki_t>
inline void operator()(problem_t &f, value_t &tn, state_t &un, array_ki_t &G, value_t &dt, state_t &unp1)#

iteration of ROCK2 method

Template Parameters:

  • problem_t – type of \(f\)

  • state_t – type of current state

  • array_ki_t – type of temporary stages (only 3 needed for ROCK2)

Parameters:

  • f – operator \(f\)

  • tn – current time

  • un – current state

  • G – array of temporary stages

  • dt – current time step

  • unp1 – returns solution at time \(t^{n+1} = t^n + \Delta t\)

template<typename rock_t = rock_coeff>
inline auto &abs_tol(value_t tol_)#

set absolute tolerance in chained config

Parameters:

tol_ – tolerance

Returns:

auto& returns this object

template<typename rock_t = rock_coeff>
inline auto &rel_tol(value_t tol_)#

set relative tolerance in chained config

Parameters:

tol_ – tolerance

Returns:

auto& returns this object

Helper functions#

template<bool is_embedded = false, typename value_t = double, typename eig_computer_t>
auto ponio::runge_kutta::rock::rock2(eig_computer_t &&eig_computer)#

helper to build a rock2_impl object

Template Parameters:

  • is_embedded – define if method is used as adaptive or constant time step method [default is false]

  • value_t – type of coefficients

  • eig_computer_t – type of computer of maximal eigenvalue

Parameters:

eig_computer – computer of maximal eigenvalues

template<bool is_embedded = false, typename value_t = double>
auto ponio::runge_kutta::rock::rock2()#

helper to build a rock2_impl object with power method to compute maximal eigen value

Template Parameters:

  • is_embedded – define if method is used as adaptive or constant time step method [default is false]

  • value_t – type of coefficients

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.

template<typename eig_computer_t, bool _is_embedded = false, typename _value_t = double>
class rock4_impl#

define ROCK4 method

Template Parameters:

  • eig_computer_t – type of computer of maximal eigenvalue

  • _is_embedded – define if method is used as adaptive or constant time step method [default is false]

  • _value_t – type of coefficients

Public Functions

inline rock4_impl(eig_computer_t &&_eig_computer)#

Construct a new ROCK4 algorithm object.

Parameters:

_eig_computer – estimator of spectral radius

template<typename state_t>
inline auto error(state_t const &unp1, state_t const &tmp)#

computes an error for adaptive time step method

Template Parameters:

state_t – type of current state

Parameters:

  • unp1 – estimation of solution at time \(t^{n+1} = t^n+\Delta t\)

  • tmp – other estimation of solution at time \(t^{n+1} = t^n+\Delta t\)

Returns:

auto estimation of error to compare to 1

template<typename problem_t, typename state_t, typename array_ki_t>
inline void operator()(problem_t &f, value_t &tn, state_t &un, array_ki_t &G, value_t &dt, state_t &unp1)#

iteration of ROCK4 method

Template Parameters:

  • problem_t – type of \(f\)

  • state_t – type of current state

  • array_ki_t – type of temporary stages (only 6 needed for ROCK4)

Parameters:

  • f – operator \(f\)

  • tn – current time

  • un – current state

  • G – array of temporary stages

  • dt – current time step

  • unp1 – solution \(u^{n+1}\) at time \(t^{n+1} = t^n + \Delta t\)

inline auto &info()#

gets iteration_info object

inline auto const &info() const#

gets iteration_info object (constant version)

template<typename rock_t = rock_coeff>
inline auto &abs_tol(value_t tol_)#

set absolute tolerance in chained config

Parameters:

tol_ – tolerance

Returns:

auto& returns this object

template<typename rock_t = rock_coeff>
inline auto &rel_tol(value_t tol_)#

set relative tolerance in chained config

Parameters:

tol_ – tolerance

Returns:

auto& returns this object

Helper functions#

template<bool is_embedded = false, typename value_t = double, typename eig_computer_t>
auto ponio::runge_kutta::rock::rock4(eig_computer_t &&eig_computer)#

helper to build a rock4_impl object

Template Parameters:

  • is_embedded – define if method is used as adaptive or constant time step method [default is false]

  • value_t – type of coefficients

  • eig_computer_t – type of computer of maximal eigenvalue

Parameters:

eig_computer – computer of maximal eigenvalues

template<bool is_embedded = false, typename value_t = double>
auto ponio::runge_kutta::rock::rock4()#

helper to build a rock4_impl object with power method to compute maximal eigen value

Template Parameters:

  • is_embedded – define if method is used as adaptive or constant time step method [default is false]

  • value_t – type of coefficients

Runge-Kutta Legendre method#

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

Runge-Kutta Legend first order method#

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.

template<std::size_t N_stages_, typename _value_t = double>
class explicit_rkl1#

define RKL1 (Runge-Kutta-Legendre method of ordre 1) with user defined number of stages

Warning

the method is only presented with an autonomous problem, ie \(\dot{y} = f(y)\)

Template Parameters:

  • N_stages_ – number of stages

  • _value_t – type of coefficients

Public Functions

template<typename problem_t, typename state_t, typename array_ki_t, std::size_t j>
inline void stage(Stage<j>, problem_t &f, value_t tn, state_t&, array_ki_t const &Yj, value_t dt, state_t &ui, state_t &yi)#

compute stage \(j\) of RKL1 method

\(y^{(j)} = \mu_j y^{(j-1)} + \nu_j y^{(j-2)} + \tilde{\mu}_j \Delta t f(t^n, y^{(j-1)})\)

Template Parameters:

  • problem_t – type of operator \(f\)

  • state_t – type of current state

  • array_ki_t – type of array of temporary stages

  • j – integer of stage

Parameters:

  • f – operator \(f\)

  • tn – current time

  • Yj – array of temporary stages

  • dt – current time step

  • ui – temporary step

  • yi – computed output step

template<typename problem_t, typename state_t, typename array_ki_t>
inline void stage(Stage<0>, problem_t&, value_t, state_t &yn, array_ki_t const&, value_t, state_t&, state_t &yi)#

compute stage \(0\) of RKL1 method

\(y^{(0)} = y^n\)

Template Parameters:

  • problem_t – type of operator \(f\)

  • state_t – type of current state

  • array_ki_t – type of array of temporary stages

Parameters:

  • yn – current state

  • yi – computed output step

template<typename problem_t, typename state_t, typename array_ki_t>
inline void stage(Stage<1>, problem_t &f, value_t tn, state_t &yn, array_ki_t const&, value_t dt, state_t &ui, state_t &yi)#

compute stage \(1\) of RKL1 method

\(y^{(1)} = y^n + \tilde{\mu}_1 \Delta t f(t^n, y^n)\)

Template Parameters:

  • problem_t – type of operator \(f\)

  • state_t – type of current state

  • array_ki_t – type of array of temporary stages

Parameters:

  • f – operator \(f\)

  • tn – current time

  • yn – current state

  • dt – current time step

  • ui – temporary step

  • yi – computed output step

inline auto &info()#

gets iteration_info object

inline auto const &info() const#

gets iteration_info object (constant version)

Public Static Functions

template<std::size_t j>
static inline value_t mu()#

compute \(\mu_j = \frac{2j-1}{j}\)

Template Parameters:

j – index \(j\)

template<std::size_t j>
static inline value_t nu()#

compute \(\nu_j = \frac{1-j}{j}\)

Template Parameters:

j – index \(j\)

template<std::size_t j>
static inline value_t mu_t()#

compute \(\tilde{\mu}_j = \frac{2j-1}{j}\frac{2}{s^2 + s}\) with \(s\) the number of stages

Template Parameters:

j – index \(j\)

Warning

In ponio library the number of stages of RKL1 is static, you can’t get dynamic number of stages for this method, look at ponio::runge_kutta::rock::rock2() or ponio::runge_kutta::rock::rock4() methods for dynamic number of stages (and optional adaptive time step method).

Runge-Kutta Legend second order method#

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

template<std::size_t N_stages_, typename _value_t = double>
class explicit_rkl2#

define RKL2 (Runge-Kutta-Legendre method of ordre 2) with user defined number of stages

Warning

the method is only presented with an autonomous problem, ie \(\dot{y} = f(y)\)

Template Parameters:

  • N_stages_ – number of stages

  • _value_t – type of coefficients

Public Functions

template<typename problem_t, typename state_t, typename array_ki_t, std::size_t j>
inline void stage(Stage<j>, problem_t &f, value_t tn, state_t &yn, array_ki_t const &Yj, value_t dt, state_t &ui, state_t &yi)#

compute stage \(j\) of RKL2 method

\(y^{(j)} = \mu_j y^{(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)}) + \gamma_j\Delta t f(t^n, y^{(0)})\)

Template Parameters:

  • problem_t – type of operator \(f\)

  • state_t – type of current state

  • array_ki_t – type of array of temporary stages

  • j – integer of stage

Parameters:

  • f – operator \(f\)

  • tn – current time

  • yn – current state

  • Yj – array of temporary stages

  • dt – current time step

  • ui – temporary step

  • yi – computed output step

template<typename problem_t, typename state_t, typename array_ki_t>
inline void stage(Stage<0>, problem_t &f, value_t tn, state_t &yn, array_ki_t const&, value_t dt, state_t &ui, state_t &yi)#

compute \(\Delta t f(t^n, y^n)\) term (needs for all other stages)

Warning

first stage of RKL2 method is \(y^n\) but here we compute \(\Delta t f(t^n, y^n)\) term

Template Parameters:

  • problem_t – type of operator \(f\)

  • state_t – type of current state

  • array_ki_t – type of array of temporary stages

Parameters:

  • f – operator \(f\)

  • tn – current time

  • yn – current state

  • dt – current time step

  • ui – temporary step

  • yi – computed output step

template<typename problem_t, typename state_t, typename array_ki_t>
inline void stage(Stage<1>, problem_t&, value_t, state_t &yn, array_ki_t const &Yj, value_t, state_t&, state_t &yi)#

compute stage \(1\) of RKL2 method

\(y^{(1)} = y^n + \tilde{\mu}_1 \Delta t f(t^n, y^n)\)

Template Parameters:

  • problem_t – type of operator \(f\)

  • state_t – type of current state

  • array_ki_t – type of array of temporary stages

Parameters:

  • yn – current state

  • Yj – array of temporary stages

  • yi – computed output step

template<typename problem_t, typename state_t, typename array_ki_t>
inline void stage(Stage<2>, problem_t &f, value_t tn, state_t &yn, array_ki_t const &Yj, value_t dt, state_t &ui, state_t &yi)#

compute stage \(2\) of RKL2 method

\(y^{(2)} = \mu_2 y^{(1)} + \nu_2y^n + (1-\mu_2-\nu_2)y^n + \tilde{\mu}_2\Delta tf(t^n, y^{(1)}) + \gamma_2\Delta t f(t^n, y^n)\)

Template Parameters:

  • problem_t – type of operator \(f\)

  • state_t – type of current state

  • array_ki_t – type of array of temporary stages

Parameters:

  • f – operator \(f\)

  • tn – current time

  • yn – current state

  • Yj – array of temporary stages

  • dt – current time step

  • ui – temporary step

  • yi – computed output step

inline auto &info()#

gets iteration_info object

inline auto const &info() const#

gets iteration_info object (constant version)

Public Static Functions

static inline value_t w1()#

compute \(w_1\) coefficient

\(w_1 = \frac{4}{s^2 + s -2}\) with \(s\) the number of stages of the method

template<std::size_t j>
static inline value_t mu()#

compute \(\mu_j\) coefficient

\(b_1 = \frac{1}{3}\)

Template Parameters:

j – index \(j\)

template<std::size_t j>
static inline value_t nu()#

compute \(\mu_j\) coefficient

\(mu_j = \frac{2j-1}{j}\frac{b_j}{b_{j-1}}\)

Template Parameters:

j – index \(j\)

template<std::size_t j>
static inline value_t mu_t()#

compute \(\tilde{\mu}_j\) coefficient

\(\tilde{\mu}_j = \mu_j w_1\) for \(1<j\), \(\tilde{\mu}_1 - b_1w_1\)

Template Parameters:

j – index \(j\)

template<std::size_t j>
static inline value_t gamma_t()#

compute \(\gamma_j\) coefficient

\(gamma_j = -a_{j-1}\tilde{\mu}_j\)

Template Parameters:

j – index \(j\)

Warning

In ponio library the number of stages of RKL2 is static, you can’t get dynamic number of stages for this method, look at ponio::runge_kutta::rock::rock2() or ponio::runge_kutta::rock::rock4() methods for dynamic number of stages (and optional adaptive time step method).