List of 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\).

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)\]
template<typename _value_t, typename ...methods_t>
class lie : public ponio::splitting::detail::splitting_base<_value_t, methods_t...>#

Lie splitting method.

Template Parameters:

  • value_t – type of time steps

  • methods_t – list of methods to solve each sub-problem

Public Functions

template<std::size_t I = 0, typename Problem_t, typename state_t>
inline void _call_inc(Problem_t &f, value_t tn, state_t &ui, value_t dt, state_t &uip1)#

incremental call of each method of each subproblem

Parameters:

  • f – problem to solve

  • tn – current time \(\Delta t\)

  • ui – initial solution for step I

  • dt – time step \(\Delta t\)

  • uip1 – solution \(\texttt{uip1} = \phi^{[i]}(\texttt{ui})\)

template<typename Problem_t, typename state_t>
void operator()(Problem_t &f, value_t &tn, state_t &un, value_t &dt, state_t &unp1)#

call operator to initiate Lie splitting recursion

Parameters:

  • f – problem to solve

  • tn – current time \(t^n\)

  • un – current solution \(u^n \approx u(t^n)\)

  • dt – time step \(\Delta t\)

  • unp1 – solution 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<std::size_t I>
inline auto &stages(sub_method<I>)#

gets array of stages of the Ith method

Template Parameters:

I – index of step method

template<std::size_t I>
inline auto const &stages(sub_method<I>) const#

gets array of stages of the Ith method (constant version)

Template Parameters:

I – index of step method

splitting_base(std::tuple<methods_t...> const &meths, std::array<value_t, sizeof...(methods_t)> const &dts)#

Construct a new splitting base<value t, methods t…>::splitting base object.

Template Parameters:

  • value_t – type of coefficients and time steps

  • methods_t – list of types of methods to solve each sub-problem

Parameters:

  • meths – list of methods to solve each sub-problem

  • dts – list of time step to solve each sub-problem (to iterate on each sub-step)

Helper function#

template<typename value_t, typename ...algorithms_t>
auto ponio::splitting::lie::make_lie_tuple(std::pair<algorithms_t, value_t>&&... args)#

a helper factory for ponio::splitting::detail::splitting_tuple from a tuple of algorithms to build a Lie method

Template Parameters:

  • value_t – type of coefficients

  • algorithms_t – variadic list of types of algorithms

Parameters:

args – variadic list of pairs of algorithm and time step

Returns:

a ponio::splitting::detail::splitting_tuple object build from the tuple of methods

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)\]
template<typename _value_t, typename ...methods_t>
class strang : public ponio::splitting::detail::splitting_base<_value_t, methods_t...>#

Strang splitting method.

Template Parameters:

methods_t – list of methods to solve each sub-problem

Public Functions

template<std::size_t I = 0, typename Problem_t, typename state_t>
inline void _call_inc(Problem_t &f, value_t tn, state_t &ui, value_t dt, state_t &uip1)#

incremental call of each method of each subproblem

Template Parameters:

I – solving step

Parameters:

  • f – problem to solve

  • tn – current time \(\Delta t\)

  • ui – initial solution for step I

  • dt – time step \(\Delta t\)

  • uip1 – solution \(\texttt{uip1} = \phi^{[i]}(\texttt{ui})\)

template<std::size_t I = N_methods - 1, typename Problem_t, typename state_t>
inline void _call_dec(Problem_t &f, value_t tn, state_t &ui, value_t dt, state_t &uip1)#

decremental call of each method of each subproblem

Template Parameters:

I – solving step

Parameters:

  • f – problem to solve

  • tn – current time \(\Delta t\)

  • ui – initial solution for step I

  • dt – time step \(\Delta t\)

  • uip1 – solution \(\texttt{uip1} = \phi^{[i]}(\texttt{ui})\)

template<typename Problem_t, typename state_t>
void operator()(Problem_t &f, value_t &tn, state_t &un, value_t &dt, state_t &unp1)#

call operator to initiate Strang splitting recursion

Parameters:

  • f – problem to solve

  • tn – current time \(t^n\)

  • un – current solution \(u^n \approx u(t^n)\)

  • dt – time step \(\Delta t\)

  • unp1 – solution 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<std::size_t I>
inline auto &stages(sub_method<I>)#

gets array of stages of the Ith method

Template Parameters:

I – index of step method

template<std::size_t I>
inline auto const &stages(sub_method<I>) const#

gets array of stages of the Ith method (constant version)

Template Parameters:

I – index of step method

splitting_base(std::tuple<methods_t...> const &meths, std::array<value_t, sizeof...(methods_t)> const &dts)#

Construct a new splitting base<value t, methods t…>::splitting base object.

Template Parameters:

  • value_t – type of coefficients and time steps

  • methods_t – list of types of methods to solve each sub-problem

Parameters:

  • meths – list of methods to solve each sub-problem

  • dts – list of time step to solve each sub-problem (to iterate on each sub-step)

Helper function#

template<typename value_t, typename ...Algorithms_t>
auto ponio::splitting::strang::make_strang_tuple(std::pair<Algorithms_t, value_t>&&... args)#

a helper factory for ponio::splitting::detail::splitting_tuple from a tuple of algorithms to build a Strang method

Template Parameters:

  • value_t – type of coefficients

  • Algorithms_t – variadic list of types of algorithms

Parameters:

args – variadic list of pairs of algorithm and time step

Returns:

a ponio::splitting::detail::splitting_tuple object build from the tuple of methods

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 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}\]
template<typename _value_t, typename ...methods_t>
class adaptive_strang : public ponio::splitting::detail::splitting_base<_value_t, methods_t...>#

adaptive time step Strang splitting method

Template Parameters:

methods_t – list of methods to solve each sub-problem

Public Functions

template<typename Problem_t, typename state_t>
inline void operator()(Problem_t &f, value_t &tn, state_t &un, value_t &dt, state_t &u_np1)#

call operator to initiate adaptive Strang splitting recursion

Parameters:

  • f – problem to solve

  • tn – current time \(t^n\)

  • un – current solution \(u^n \approx u(t^n)\)

  • dt – time step \(\Delta t\)

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

template<typename Problem_t, typename state_t>
inline auto lipschitz_constant_estimate(Problem_t &f, value_t tn, state_t &un, value_t dt)#

estimate Lipschitz constant of problem \(f\) at time \(t^n\)

The error of classical splitting method is:

\[ \mathcal{S}^{\Delta t}u_0 - \mathcal{T}^{\Delta t}u_0 = C_0 \Delta t^3 \]

where \(\mathcal{S}\) is the time flow of Strang method, and \(\mathcal{T}\) is the exact time flow of the equation, with a shifted splitting method, we can compute

\[ \mathcal{S}^{\Delta t}u_0 - \mathcal{S}_\delta^{\Delta t} u_0 = \delta C_\delta \Delta t^2 \]

we find the optimal time step \(\Delta t^\star\) as

\[ \| \mathcal{S}^{\Delta t}u_0 - \mathcal{T}^{\Delta t}u_0 \| \leq \| \mathcal{S}^{\Delta t}u_0 - \mathcal{S}_\delta^{\Delta t}u_0 \|,\ \forall \Delta t \leq \Delta t^\star \]

we get an estimate of \(\Delta t^\star\)

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

To estimate \(C_0\) we need compute two errors

\[\begin{split} \begin{aligned} e_1 &= \| \mathcal{S}^{a_1\Delta t}u_0 \mathcal{S}^{b_1\Delta t}\mathcal{S}^{c_1\Delta t}u_0 \| \\ e_2 &= \| \mathcal{S}^{a_2\Delta t}u_0 \mathcal{S}^{b_2\Delta t}\mathcal{S}^{c_2\Delta t}u_0 \| \end{aligned} \end{split}\]

we obtain two inequalities

\[\begin{split} \begin{aligned} \|e_1 - (a_1^3 - b_1^3)C_0\Delta t^3\| \leq \omega C_0 c_1^3\Delta t^3 \\ \|e_2 - (a_2^3 - b_2^3)C_0\Delta t^3\| \leq \omega C_0 c_2^3\Delta t^3 \end{aligned} \end{split}\]

we choose parameters as

\[ a_1 = 1,\quad b_1=\frac{1}{2},\quad c_1=\frac{1}{2},\quad a_2=c_1,\quad b_2=\frac{2}{5},\quad c_2=\frac{1}{10} \]

to obtain estimate of \(C_0\) (error estimate) and \(\omega\) (Lipschitz constant estimate).

Parameters:

  • f – callable obect which represents the problem to solve

  • tn – time \(t^n\) last time where solution is computed

  • un – computed solution \(u^n\) à time \(t^n\)

  • dt – time step

inline auto &info()#

gets iteration_info object

inline auto const &info() const#

gets iteration_info object (constant version)

template<std::size_t I>
inline auto &stages(sub_method<I>)#

gets array of stages of the Ith method

Template Parameters:

I – index of step method

template<std::size_t I>
inline auto const &stages(sub_method<I>) const#

gets array of stages of the Ith method (constant version)

Template Parameters:

I – index of step method

splitting_base(std::tuple<methods_t...> const &meths, std::array<value_t, sizeof...(methods_t)> const &dts)#

Construct a new splitting base<value t, methods t…>::splitting base object.

Template Parameters:

  • value_t – type of coefficients and time steps

  • methods_t – list of types of methods to solve each sub-problem

Parameters:

  • meths – list of methods to solve each sub-problem

  • dts – list of time step to solve each sub-problem (to iterate on each sub-step)

Helper function#

template<typename value_t, typename ...Algorithms_t>
auto ponio::splitting::strang::make_adaptive_strang_tuple(value_t delta, value_t tol, std::pair<Algorithms_t, value_t>&&... args)#

a helper factory for ponio::splitting::detail::splitting_tuple from a tuple of algorithms to build an adaptive time step Strang method

Template Parameters:

  • value_t – type of coefficients

  • Algorithms_t – variadic list of types of algorithms

Parameters:

  • delta – shift argument

  • tol – tolerance for adaptive time step algorithm

  • args – variadic list of pairs of algorithm and time step

Returns:

a ponio::splitting::detail::splitting_tuple object build from the tuple of methods

Generic splitting tuple#

The class ponio::splitting::detail::splitting_tuple is useful to build each splitting method from variadic number of pair of algorithm and time step. This class can has optional arguments, stored in a tuple, to call constructor of _splitting_method_t (which is ponio::splitting::lie::lie, ponio::splitting::strang::strang or ponio::splitting::strang::adaptive_strang).

template<template<typename, typename...> typename _splitting_method_t, typename value_t, typename optional_args_t, typename ...Algorithms_t>
class splitting_tuple#

generic class to prevent code duplication between splitting method

Template Parameters:

  • _splitting_method_t – type of splitting method

  • value_t – type of coefficient

  • optional_args_t – type of tuple of optional arguments (void if not needed)

  • Algorithms_t – type of algorithms to solve each step of splitting

Helper functions#

template<template<typename, typename...> typename _splitting_method_t, typename value_t, typename ...Methods_t>
auto ponio::splitting::detail::make_splitting_from_tuple(std::tuple<Methods_t...> const &meths, std::array<value_t, sizeof...(Methods_t)> const &dts)#

factory for generic splitting method (strang, lie)

Template Parameters:

  • _splitting_method_t – type of splitting method

  • value_t – type of coefficients

  • Methods_t – type of methods to solve each step of splitting

Parameters:

  • meths – tuple of methods

  • dts – time step for each method

template<template<typename, typename...> typename _splitting_method_t, typename value_t, typename optional_tuple_t, typename ...Methods_t>
auto ponio::splitting::detail::make_splitting_from_tuple(std::tuple<Methods_t...> const &meths, std::array<value_t, sizeof...(Methods_t)> const &dts, optional_tuple_t optional_args)#

factory for generic splitting method (adaptive_strang)

Template Parameters:

  • _splitting_method_t – type of splitting method

  • value_t – type of coefficients

  • optional_tuple_t – type of tuple of optional arguments

  • Methods_t – type of methods to solve each step of splitting

Parameters:

  • meths – tuple of methods

  • dts – time step for each method

  • optional_args – tuple of optional arguments