ODE integration

\( \dot y = f(y,t) \) More...

Collaboration diagram for ODE integration:

Topics
	Utilities for ODE integrators

Classes
struct	dg::Adaptive< Stepper >
	Driver class for adaptive timestep ODE integration. More...
struct	dg::AdaptiveTimeloop< ContainerType >
	Integrate using a while loop. More...
struct	dg::ExplicitMultistep< ContainerType >
	General explicit linear multistep ODE integrator \( \begin{align} v^{n+1} = \sum_{j=0}^{s-1} a_j v^{n-j} + \Delta t\left(\sum_{j=0}^{s-1}b_j \hat f\left(t^{n}-j\Delta t, v^{n-j}\right)\right) \end{align} \). More...
struct	dg::ImExMultistep< ContainerType >
	Semi-implicit multistep ODE integrator \( \begin{align} v^{n+1} = \sum_{q=0}^{s-1} a_q v^{n-q} + \Delta t\left[\left(\sum_{q=0}^{s-1}b_q \hat E(t^{n}-q\Delta t, v^{n-q}) + \sum_{q=1}^{s} c_q \hat I( t^n - q\Delta t, v^{n-q})\right) + c_0\hat I(t^{n}+\Delta t, v^{n+1})\right] \end{align} \). More...
struct	dg::ImplicitMultistep< ContainerType >
	Implicit multistep ODE integrator \( \begin{align} v^{n+1} &= \sum_{i=0}^{s-1} a_i v^{n-i} + \Delta t \sum_{i=1}^{s} c_i\hat I(t^{n+1-i}, v^{n+1-i}) + \Delta t c_{0} \hat I (t + \Delta t, v^{n+1}) \\ \end{align} \). More...
struct	dg::FilteredExplicitMultistep< ContainerType >
	EXPERIMENTAL: General explicit linear multistep ODE integrator with Limiter / Filter \( \begin{align} \tilde v &= \sum_{j=0}^{s-1} a_j v^{n-j} + \Delta t\left(\sum_{j=0}^{s-1}b_j \hat f\left(t^{n}-j\Delta t, v^{n-j}\right)\right) \\ v^{n+1} &= \Lambda\Pi \left( \tilde v\right) \end{align} \). More...
struct	dg::MultistepTimeloop< ContainerType >
	Integrate using a for loop and a fixed non-changeable time-step. More...
struct	dg::aTimeloop< ContainerType >
	Abstract timeloop independent of stepper and ODE. More...
struct	dg::ERKStep< ContainerType >
	Embedded Runge Kutta explicit time-step with error estimate \( \begin{align} k_i = f\left( t^n + c_i \Delta t, u^n + \Delta t \sum_{j=1}^{i-1} a_{ij} k_j\right) \\ u^{n+1} = u^{n} + \Delta t\sum_{j=1}^s b_j k_j \\ \tilde u^{n+1} = u^{n} + \Delta t\sum_{j=1}^s \tilde b_j k_j \\ \delta^{n+1} = u^{n+1} - \tilde u^{n+1} = \Delta t\sum_{j=1}^s (b_j - \tilde b_j) k_j \end{align} \). More...
struct	dg::FilteredERKStep< ContainerType >
	EXPERIMENTAL: Filtered Embedded Runge Kutta explicit time-step with error estimate \( \begin{align} k_i = f\left( t^n + c_i \Delta t, \Lambda\Pi \left[u^n + \Delta t \sum_{j=1}^{i-1} a_{ij} k_j\right]\right) \\ u^{n+1} = \Lambda\Pi\left[u^{n} + \Delta t\sum_{j=1}^s b_j k_j\right] \\ \delta^{n+1} = \Delta t\sum_{j=1}^s (\tilde b_j - b_j) k_j \end{align} \). More...
struct	dg::ARKStep< ContainerType >
	Additive Runge Kutta (semi-implicit) time-step with error estimate following The ARKode library More...
struct	dg::DIRKStep< ContainerType >
	Embedded diagonally implicit Runge Kutta time-step with error estimate \( \begin{align} k_i = f\left( t^n + c_i \Delta t, u^n + \Delta t \sum_{j=1}^{i} a_{ij} k_j\right) \\ u^{n+1} = u^{n} + \Delta t\sum_{j=1}^s b_j k_j \\ \tilde u^{n+1} = u^{n} + \Delta t\sum_{j=1}^s \tilde b_j k_j \\ \delta^{n+1} = u^{n+1} - \tilde u^{n+1} = \Delta t\sum_{j=1}^s (b_j-\tilde b_j) k_j \end{align} \). More...
struct	dg::ShuOsher< ContainerType >
	Shu-Osher fixed-step explicit ODE integrator with Slope Limiter / Filter \( \begin{align} u_0 &= u_n \\ u_i &= \Lambda\Pi \left(\sum_{j=0}^{i-1}\left[ \alpha_{ij} u_j + \Delta t \beta_{ij} f( t_j, u_j)\right]\right)\\ u^{n+1} &= u_s \end{align} \). More...
struct	dg::SinglestepTimeloop< ContainerType >
	Integrate using a for loop and a fixed time-step. More...

Typedefs
template<class ContainerType >
using	dg::RungeKutta = ERKStep<ContainerType>
	Runge-Kutta fixed-step explicit ODE integrator \( \begin{align} k_i = f\left( t^n + c_i \Delta t, u^n + \Delta t \sum_{j=1}^{i-1} a_{ij} k_j\right) \\ u^{n+1} = u^{n} + \Delta t\sum_{j=1}^s b_j k_j \end{align} \).
template<class ContainerType >
using	dg::FilteredRungeKutta = FilteredERKStep<ContainerType>
	Filtered Runge-Kutta fixed-step explicit ODE integrator \( \begin{align} k_i = f\left( t^n + c_i \Delta t, \Lambda\Pi \left[u^n + \Delta t \sum_{j=1}^{i-1} a_{ij} k_j\right]\right) \\ u^{n+1} = \Lambda\Pi\left[u^{n} + \Delta t\sum_{j=1}^s b_j k_j\right] \end{align} \).
template<class ContainerType >
using	dg::ImplicitRungeKutta = DIRKStep<ContainerType>
	Runge-Kutta fixed-step implicit ODE integrator \( \begin{align} k_i = f\left( t^n + c_i \Delta t, u^n + \Delta t \sum_{j=1}^{s} a_{ij} k_j\right) \\ u^{n+1} = u^{n} + \Delta t\sum_{j=1}^s b_j k_j \end{align} \).

Detailed Description

\( \dot y = f(y,t) \)

Typedef Documentation

FilteredRungeKutta

template<class ContainerType >

using dg::FilteredRungeKutta = FilteredERKStep<ContainerType>

Filtered Runge-Kutta fixed-step explicit ODE integrator \( \begin{align} k_i = f\left( t^n + c_i \Delta t, \Lambda\Pi \left[u^n + \Delta t \sum_{j=1}^{i-1} a_{ij} k_j\right]\right) \\ u^{n+1} = \Lambda\Pi\left[u^{n} + \Delta t\sum_{j=1}^s b_j k_j\right] \end{align} \).

The method is defined by its (explicit) ButcherTableau, given by the coefficients a, b and c, and s is the number of stages.

You can provide your own coefficients or use one of our predefined methods (including the ones in Shu-Osher form):

We follow the naming convention of the ARKode library https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html (They also provide nice stability plots for their methods) as NAME-S-P-Q or NAME-S-Q, where

• NAME is the author or name of the method
• S is the number of stages in the method
• P is the global order of the embedding
• Q is the global order of the method

Name	Identifier	Description
Euler	dg::EXPLICIT_EULER_1_1	Explicit Euler
Midpoint-2-2	dg::MIDPOINT_2_2	Midpoint method
Kutta-3-3	dg::KUTTA_3_3	Kutta-3-3
Runge-Kutta-4-4	dg::CLASSIC_4_4	"The" Runge-Kutta method
SSPRK-2-2	dg::SSPRK_2_2	SSPRK (2,2) CFL_eff = 0.5 "Shu-Osher-Form"
SSPRK-3-2	dg::SSPRK_3_2	SSPRK (3,2) CFL_eff = 0.66 "Shu-Osher-Form"
SSPRK-3-3	dg::SSPRK_3_3	SSPRK (3,3) CFL_eff = 0.33 "Shu-Osher-Form"
SSPRK-5-3	dg::SSPRK_5_3	SSPRK (5,3) CFL_eff = 0.5 "Shu-Osher-Form"
SSPRK-5-4	dg::SSPRK_5_4	SSPRK (5,4) CFL_eff = 0.37 "Shu-Osher-Form"
Heun-Euler-2-1-2	dg::HEUN_EULER_2_1_2	Heun-Euler-2-1-2
Cavaglieri-3-1-2 (explicit)	dg::CAVAGLIERI_3_1_2	Low-storage implicit/explicit Runge-Kutta schemes for the simulation of stiff high-dimensional ODE systems IMEXRKCB2 scheme
Fehlberg-3-2-3	dg::FEHLBERG_3_2_3	The original uses the embedding as the solution (but we do not) [Hairer, Noersett, Wanner, Solving ordinary differential Equations I, 1987]
Bogacki-Shampine-4-2-3	dg::BOGACKI_SHAMPINE_4_2_3	Bogacki-Shampine (fsal)
Cavaglieri-4-2-3 (explicit)	dg::CAVAGLIERI_4_2_3	Low-storage implicit/explicit Runge-Kutta schemes for the simulation of stiff high-dimensional ODE systems The SSP scheme IMEXRKCB3c (explicit part)
ARK-4-2-3 (explicit)	dg::ARK324L2SA_ERK_4_2_3	ARK-4-2-3 (explicit)
Zonneveld-5-3-4	dg::ZONNEVELD_5_3_4	Zonneveld-5-3-4
ARK-6-3-4 (explicit)	dg::ARK436L2SA_ERK_6_3_4	ARK-6-3-4 (explicit)
Sayfy_Aburub-6-3-4	dg::SAYFY_ABURUB_6_3_4	Sayfy_Aburub_6_3_4
Cash_Karp-6-4-5	dg::CASH_KARP_6_4_5	Cash-Karp
Fehlberg-6-4-5	dg::FEHLBERG_6_4_5	The original uses the embedding as the solution (but we do not) Runge-Kutta-Fehlberg
Dormand-Prince-7-4-5	dg::DORMAND_PRINCE_7_4_5	Dormand-Prince method (fsal)
Tsitouras09-7-4-5	dg::TSITOURAS09_7_4_5	Tsitouras 5(4) method from 2009 (fsal), The default method in Julia
Tsitouras11-7-4-5	dg::TSITOURAS11_7_4_5	Tsitouras 5(4) method from 2011 (fsal) Further improves Tsitouras09 (Note that in the paper b-bt is printed instead of bt)
ARK-8-4-5 (explicit)	dg::ARK548L2SA_ERK_8_4_5	Kennedy and Carpenter (2019) Optimum ARK_2 method (explicit part)
Verner-9-5-6	dg::VERNER_9_5_6	Verner-9-5-6 (fsal)
Verner-10-6-7	dg::VERNER_10_6_7	Verner-10-6-7
Fehlberg-13-7-8	dg::FEHLBERG_13_7_8	Fehlberg-13-7-8
Dormand-Prince-13-7-8	dg::DORMAND_PRINCE_13_7_8	[Hairer, Noersett, Wanner, Solving ordinary differential Equations I, 1987]
Feagin-17-8-10	dg::FEAGIN_17_8_10	Feagin (The RK10(8) method)

The following code snippet demonstrates how to use the class for the integration of the harmonic oscillator:

    auto rhs = [=]( double t, const std::array<double,2>& y,
            std::array<double,2>& yp)
    {
        //damped driven harmonic oscillator
        // x -> y[0] , v -> y[1]
        yp[0] = y[1];
        yp[1] = -2.*damping*omega_0*y[1] - omega_0*omega_0*y[0]
                + sin(omega_drive*t);
    };
    auto solve = [=]( double alpha, double t, std::array<double,2>& y,
            const std::array<double,2>& yp)
    {
        // y - alpha RHS( t, y) = rho
        // can be solved analytically
        y[1] = ( yp[1] + alpha*sin(omega_drive*t) - alpha*omega_0*omega_0*yp[0])/
               (1.+2.*alpha*damping*omega_0+alpha*alpha*omega_0*omega_0);
        y[0] = yp[0] + alpha*y[1];
    };
    SECTION( "Runge-Kutta-4-4")
    {

    //... in main
    //set start and end time, number of steps and timestep
    const double t_start = 0., t_end = 1.;
    const unsigned N = 20;
    const double dt = (t_end - t_start)/(double)N;
    //set physical parameters and initial condition
    const double damping = 0.2, omega_0 = 1.0, omega_drive = 0.9;
    std::array<double,2> u = solution(t_start, damping, omega_0, omega_drive);
    //construct a functor with the right interface
    auto rhs = [=]( double t, const std::array<double,2>& y,
            std::array<double,2>& yp)
    {
        //damped driven harmonic oscillator
        // x -> y[0] , v -> y[1]
        yp[0] = y[1];
        yp[1] = -2.*damping*omega_0*y[1] - omega_0*omega_0*y[0]
                + sin(omega_drive*t);
    };
    auto solve = [=]( double alpha, double t, std::array<double,2>& y,
            const std::array<double,2>& yp)
    {
        // y - alpha RHS( t, y) = rho
        // can be solved analytically
        y[1] = ( yp[1] + alpha*sin(omega_drive*t) - alpha*omega_0*omega_0*yp[0])/
               (1.+2.*alpha*damping*omega_0+alpha*alpha*omega_0*omega_0);
        y[0] = yp[0] + alpha*y[1];
    };
    SECTION( "Runge-Kutta-4-4")
    {
        //construct Runge Kutta class
        dg::RungeKutta<std::array<double,2> >  rk( "Runge-Kutta-4-4", u);
        //integration loop
        double t=t_start;
        for( unsigned i=0; i<N; i++)
            rk.step( rhs, t, u, t, u, dt); //step inplace
        //now compute error
        dg::blas1::axpby( 1., solution(t_end, damping, omega_0, omega_drive), -1., u);
        double err = sqrt( dg::blas1::dot( u,u));
        INFO( "Norm of error is "<<err);
        CHECK( err < 1e-7);
    }

Note

Uses only dg::blas1 routines to integrate one step.

Template Parameters

ContainerType

Any class for which a specialization of TensorTraits exists and which fulfills the requirements of the there defined data and execution policies derived from AnyVectorTag and AnyPolicyTag. Among others dg::HVec (serial), dg::DVec (cuda / omp), dg::MHVec (mpi + serial) or dg::MDVec (mpi + cuda / omp) std::vector<dg::DVec> (vector of shared device vectors), std::array<double, 4> (array of 4 doubles) or std::map < std::string, dg::DVec> ( a map of named vectors) double (scalar) and other primitive types ... If there are several ContainerTypes in the argument list, then TensorTraits must exist for all of them

See also

See The dg dispatch system for a detailed explanation of our type dispatch system

ImplicitRungeKutta

template<class ContainerType >

using dg::ImplicitRungeKutta = DIRKStep<ContainerType>

Runge-Kutta fixed-step implicit ODE integrator \( \begin{align} k_i = f\left( t^n + c_i \Delta t, u^n + \Delta t \sum_{j=1}^{s} a_{ij} k_j\right) \\ u^{n+1} = u^{n} + \Delta t\sum_{j=1}^s b_j k_j \end{align} \).

The method is defined by its (implicit) ButcherTableau, given by the coefficients a, b and c, and s is the number of stages.

You can provide your own coefficients or use one of our predefined methods:

We follow the naming convention of the ARKode library https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html (They also provide nice stability plots for their methods) as NAME-S-P-Q or NAME-S-Q, where

• NAME is the author or name of the method
• S is the number of stages in the method
• P is the global order of the embedding
• Q is the global order of the method

Name	Identifier	Description
Euler (implicit)	dg::IMPLICIT_EULER_1_1	backward Euler \(a_{11} = 1\)
Midpoint (implicit)	dg::IMPLICIT_MIDPOINT_1_2	implicit Midpoint \( a_{11} = 0.5 \)
Trapezoidal-2-2	dg::TRAPEZOIDAL_2_2	Crank-Nicolson method \( a_{11} = 0\ a_{22} = 0.5 \)
SDIRK-2-1-2	dg::SDIRK_2_1_2	generic 2nd order A and L-stable \( a_{ii} = 0.29 \)
Cavaglieri-3-1-2 (implicit)	dg::CAVAGLIERI_IMPLICIT_3_1_2	Low-storage implicit/explicit Runge-Kutta schemes for the simulation of stiff high-dimensional ODE systems IMEXRKCB2 scheme \( a_{11} = 0 \)
Sanchez-3-3	dg::SANCHEZ_3_3	symplectic DIRK
Kvaerno-4-2-3	dg::KVAERNO_4_2_3	Kvaerno-4-2-3 \( a_{11} = 0\ a_{ii} = 0.44 \)
SDIRK-4-2-3	dg::SDIRK_4_2_3	Cameron2002 \( a_{ii}=0.25\)
Cavaglieri-4-2-3 (implicit)	dg::CAVAGLIERI_IMPLICIT_4_2_3	Low-storage implicit/explicit Runge-Kutta schemes for the simulation of stiff high-dimensional ODE systems The SSP scheme IMEXRKCB3c (implicit part) \( a_{11} = 0 \)
ARK-4-2-3 (implicit)	dg::ARK324L2SA_DIRK_4_2_3	ARK-4-2-3 (implicit) \( a_{11} = 0\ a_{ii} = 0.44 \)
Sanchez-3-4	dg::SANCHEZ_3_4	symplectic DIRK
Cash-5-2-4	dg::CASH_5_2_4	Cash-5-2-4 \( a_{ii} = 0.44 \)
Cash-5-3-4	dg::CASH_5_3_4	Cash-5-3-4 \( a_{ii} = 0.44 \)
SDIRK-5-3-4	dg::SDIRK_5_3_4	SDIRK-5-3-4 \( a_{ii} = 0.25\)
Kvaerno-5-3-4	dg::KVAERNO_5_3_4	Kvaerno-5-3-4 \( a_{11} = 0\ a_{ii} = 0.44\)
ARK-6-3-4 (implicit)	dg::ARK436L2SA_DIRK_6_3_4	ARK-6-3-4 (implicit) \( a_{11} = 0\ a_{ii} = 0.25\)
Sanchez-6-5	dg::SANCHEZ_6_5	symplectic DIRK
Kvaerno-7-4-5	dg::KVAERNO_7_4_5	Kvaerno-7-4-5 \( a_{11} = 0 \ a_{ii} = 0.26\)
ARK-8-4-5 (implicit)	dg::ARK548L2SA_DIRK_8_4_5	Kennedy and Carpenter (2019) Optimum ARK_2 method (implicit part) \( a_{11} = 0\ a_{ii} = 0.22\)
Sanchez-7-6	dg::SANCHEZ_7_6	symplectic DIRK

Note

Schemes marked with \( a_{11} = 0\) call the implicit_rhs instead of a solve once per step (the first stage is explicit, typically named EDIRK), all others never call implicit_rhs. Schemes marked with \(a_{ii} = ...\) (lower is better) have all diagonal elements equal (typically named SDIRK for "singly"). Schemes marked with both have equal elements only for i>1 (named ESDIRK).

Uses only dg::blas1 routines to integrate one step.

Template Parameters

ContainerType

Any class for which a specialization of TensorTraits exists and which fulfills the requirements of the there defined data and execution policies derived from AnyVectorTag and AnyPolicyTag. Among others dg::HVec (serial), dg::DVec (cuda / omp), dg::MHVec (mpi + serial) or dg::MDVec (mpi + cuda / omp) std::vector<dg::DVec> (vector of shared device vectors), std::array<double, 4> (array of 4 doubles) or std::map < std::string, dg::DVec> ( a map of named vectors) double (scalar) and other primitive types ... If there are several ContainerTypes in the argument list, then TensorTraits must exist for all of them

See also

See The dg dispatch system for a detailed explanation of our type dispatch system

RungeKutta

template<class ContainerType >

using dg::RungeKutta = ERKStep<ContainerType>

Runge-Kutta fixed-step explicit ODE integrator \( \begin{align} k_i = f\left( t^n + c_i \Delta t, u^n + \Delta t \sum_{j=1}^{i-1} a_{ij} k_j\right) \\ u^{n+1} = u^{n} + \Delta t\sum_{j=1}^s b_j k_j \end{align} \).

The method is defined by its (explicit) ButcherTableau, given by the coefficients a, b and c, and s is the number of stages.

You can provide your own coefficients or use one of our predefined methods (including the ones in Shu-Osher form):

We follow the naming convention of the ARKode library https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html (They also provide nice stability plots for their methods) as NAME-S-P-Q or NAME-S-Q, where

• NAME is the author or name of the method
• S is the number of stages in the method
• P is the global order of the embedding
• Q is the global order of the method

Name	Identifier	Description
Euler	dg::EXPLICIT_EULER_1_1	Explicit Euler
Midpoint-2-2	dg::MIDPOINT_2_2	Midpoint method
Kutta-3-3	dg::KUTTA_3_3	Kutta-3-3
Runge-Kutta-4-4	dg::CLASSIC_4_4	"The" Runge-Kutta method
SSPRK-2-2	dg::SSPRK_2_2	SSPRK (2,2) CFL_eff = 0.5 "Shu-Osher-Form"
SSPRK-3-2	dg::SSPRK_3_2	SSPRK (3,2) CFL_eff = 0.66 "Shu-Osher-Form"
SSPRK-3-3	dg::SSPRK_3_3	SSPRK (3,3) CFL_eff = 0.33 "Shu-Osher-Form"
SSPRK-5-3	dg::SSPRK_5_3	SSPRK (5,3) CFL_eff = 0.5 "Shu-Osher-Form"
SSPRK-5-4	dg::SSPRK_5_4	SSPRK (5,4) CFL_eff = 0.37 "Shu-Osher-Form"
Heun-Euler-2-1-2	dg::HEUN_EULER_2_1_2	Heun-Euler-2-1-2
Cavaglieri-3-1-2 (explicit)	dg::CAVAGLIERI_3_1_2	Low-storage implicit/explicit Runge-Kutta schemes for the simulation of stiff high-dimensional ODE systems IMEXRKCB2 scheme
Fehlberg-3-2-3	dg::FEHLBERG_3_2_3	The original uses the embedding as the solution (but we do not) [Hairer, Noersett, Wanner, Solving ordinary differential Equations I, 1987]
Bogacki-Shampine-4-2-3	dg::BOGACKI_SHAMPINE_4_2_3	Bogacki-Shampine (fsal)
Cavaglieri-4-2-3 (explicit)	dg::CAVAGLIERI_4_2_3	Low-storage implicit/explicit Runge-Kutta schemes for the simulation of stiff high-dimensional ODE systems The SSP scheme IMEXRKCB3c (explicit part)
ARK-4-2-3 (explicit)	dg::ARK324L2SA_ERK_4_2_3	ARK-4-2-3 (explicit)
Zonneveld-5-3-4	dg::ZONNEVELD_5_3_4	Zonneveld-5-3-4
ARK-6-3-4 (explicit)	dg::ARK436L2SA_ERK_6_3_4	ARK-6-3-4 (explicit)
Sayfy_Aburub-6-3-4	dg::SAYFY_ABURUB_6_3_4	Sayfy_Aburub_6_3_4
Cash_Karp-6-4-5	dg::CASH_KARP_6_4_5	Cash-Karp
Fehlberg-6-4-5	dg::FEHLBERG_6_4_5	The original uses the embedding as the solution (but we do not) Runge-Kutta-Fehlberg
Dormand-Prince-7-4-5	dg::DORMAND_PRINCE_7_4_5	Dormand-Prince method (fsal)
Tsitouras09-7-4-5	dg::TSITOURAS09_7_4_5	Tsitouras 5(4) method from 2009 (fsal), The default method in Julia
Tsitouras11-7-4-5	dg::TSITOURAS11_7_4_5	Tsitouras 5(4) method from 2011 (fsal) Further improves Tsitouras09 (Note that in the paper b-bt is printed instead of bt)
ARK-8-4-5 (explicit)	dg::ARK548L2SA_ERK_8_4_5	Kennedy and Carpenter (2019) Optimum ARK_2 method (explicit part)
Verner-9-5-6	dg::VERNER_9_5_6	Verner-9-5-6 (fsal)
Verner-10-6-7	dg::VERNER_10_6_7	Verner-10-6-7
Fehlberg-13-7-8	dg::FEHLBERG_13_7_8	Fehlberg-13-7-8
Dormand-Prince-13-7-8	dg::DORMAND_PRINCE_13_7_8	[Hairer, Noersett, Wanner, Solving ordinary differential Equations I, 1987]
Feagin-17-8-10	dg::FEAGIN_17_8_10	Feagin (The RK10(8) method)

The following code snippet demonstrates how to use the class for the integration of the harmonic oscillator:

    auto rhs = [=]( double t, const std::array<double,2>& y,
            std::array<double,2>& yp)
    {
        //damped driven harmonic oscillator
        // x -> y[0] , v -> y[1]
        yp[0] = y[1];
        yp[1] = -2.*damping*omega_0*y[1] - omega_0*omega_0*y[0]
                + sin(omega_drive*t);
    };
    auto solve = [=]( double alpha, double t, std::array<double,2>& y,
            const std::array<double,2>& yp)
    {
        // y - alpha RHS( t, y) = rho
        // can be solved analytically
        y[1] = ( yp[1] + alpha*sin(omega_drive*t) - alpha*omega_0*omega_0*yp[0])/
               (1.+2.*alpha*damping*omega_0+alpha*alpha*omega_0*omega_0);
        y[0] = yp[0] + alpha*y[1];
    };
    SECTION( "Runge-Kutta-4-4")
    {

    //... in main
    //set start and end time, number of steps and timestep
    const double t_start = 0., t_end = 1.;
    const unsigned N = 20;
    const double dt = (t_end - t_start)/(double)N;
    //set physical parameters and initial condition
    const double damping = 0.2, omega_0 = 1.0, omega_drive = 0.9;
    std::array<double,2> u = solution(t_start, damping, omega_0, omega_drive);
    //construct a functor with the right interface
    auto rhs = [=]( double t, const std::array<double,2>& y,
            std::array<double,2>& yp)
    {
        //damped driven harmonic oscillator
        // x -> y[0] , v -> y[1]
        yp[0] = y[1];
        yp[1] = -2.*damping*omega_0*y[1] - omega_0*omega_0*y[0]
                + sin(omega_drive*t);
    };
    auto solve = [=]( double alpha, double t, std::array<double,2>& y,
            const std::array<double,2>& yp)
    {
        // y - alpha RHS( t, y) = rho
        // can be solved analytically
        y[1] = ( yp[1] + alpha*sin(omega_drive*t) - alpha*omega_0*omega_0*yp[0])/
               (1.+2.*alpha*damping*omega_0+alpha*alpha*omega_0*omega_0);
        y[0] = yp[0] + alpha*y[1];
    };
    SECTION( "Runge-Kutta-4-4")
    {
        //construct Runge Kutta class
        dg::RungeKutta<std::array<double,2> >  rk( "Runge-Kutta-4-4", u);
        //integration loop
        double t=t_start;
        for( unsigned i=0; i<N; i++)
            rk.step( rhs, t, u, t, u, dt); //step inplace
        //now compute error
        dg::blas1::axpby( 1., solution(t_end, damping, omega_0, omega_drive), -1., u);
        double err = sqrt( dg::blas1::dot( u,u));
        INFO( "Norm of error is "<<err);
        CHECK( err < 1e-7);
    }

Note

Uses only dg::blas1 routines to integrate one step.

Template Parameters

ContainerType

Any class for which a specialization of TensorTraits exists and which fulfills the requirements of the there defined data and execution policies derived from AnyVectorTag and AnyPolicyTag. Among others dg::HVec (serial), dg::DVec (cuda / omp), dg::MHVec (mpi + serial) or dg::MDVec (mpi + cuda / omp) std::vector<dg::DVec> (vector of shared device vectors), std::array<double, 4> (array of 4 doubles) or std::map < std::string, dg::DVec> ( a map of named vectors) double (scalar) and other primitive types ... If there are several ContainerTypes in the argument list, then TensorTraits must exist for all of them

See also

See The dg dispatch system for a detailed explanation of our type dispatch system