ODE Solvers¶

HighFM support several methods for solving systems of differential equations that have the following form

\[\frac{d \vec{x}}{dt} = f(t, \vec{x}), \qquad \vec{x}(t) \in \mathbb{C}^n\]

with a particular focus on systems that arise from the spatial discretisation of partial differential equations (PDEs).

Defining the equations¶

First, you need to define the system of equations using the built-in structs:

template<Number T, DataLocation L = Host> struct SemiDiscreteODE1¶

A system of 1st order differential equations. It has the following form

\[ \frac{d \vec{x}}{d t} = \alpha \mathbf{A}\vec{x} + \vec{f}(t) \]

where \( \mathbf{A} \) is a sparse matrix from a spatial discretization of a PDE and \( \vec{f}(t) \) is the forcing term.

template<Number T, DataLocation L = Host> struct SemiDiscreteODE2¶

A system of 2nd order ordinary differential equations. It has the following form

\[ \frac{d^2 \vec{x}}{d t^2} = \alpha \mathbf{A}\vec{x} + \vec{f}(t) \]

where \( \mathbf{A} \) is a sparse matrix from a spatial discretization of a PDE and \( \vec{f}(t) \) is the forcing term. Let \( \vec{y} = \frac{d}{d t} \vec{x} \), this system can be rewritten as

\[\begin{split} \frac{d}{d t} \begin{bmatrix} \vec{x} \\ \vec{y} \end{bmatrix} = \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ \alpha \mathbf{A} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \vec{x} \\ \vec{y} \end{bmatrix} + \begin{bmatrix} \vec{0} \\ \vec{f}(t) \end{bmatrix} \end{split}\]

In this form, the initial conditions must be stored as \( [\vec{x}(t_0) \; \vec{y}(t_0)]^T \), where \( t_0 \) is the initial time.

Or, alternatively, create your own struct that overloads the operator () as follows:

template<T>
struct MyODE {
    void operator()(VectorMap<T> x, value_type t, VectorMap<T> dxdt) const;
};

Here, x and t are the solution and the time of the current step, respectively. The operator () then evaluates the right-hand side f(t, x) and stores the result in dxdt.

Defining the problem¶

Next, you create a new ODEProblem object that defines the problem related to the system of equations.

template<typename System> struct ODEProblem¶

An object representing the ODE problem.

Public Members

std::array<double, 2> time_span¶: Set the initial and end time.

double time_step = kAuto ¶: Define the initial time step. If set to kAuto, the library will use a heuristic to calculate the initial time based on atol and rtol.

double max_time_step = 1.0¶: Define the maximum time step allowed.

double min_time_step = 0.0¶: Define the minimum time step allowed.

FinalTimeHandle final_time_mode = kAdaptTimeStep ¶: Determines how the last time step will be handled.

System system¶: Systems of differential equations. It can be either SemiDiscreteODE1, SemiDiscreteODE2 or a user-defined system (such that calling system() will evalute the right-hand side for a given time t).

double atol = 1E-6¶: Absolute tolerance. This parameter is not used if the ODE solver has constant step size and time_step != kAuto.

double rtol = 1E-6¶: Relative tolerance. This parameter is not used if the ODE solver has constant step size and time_step != kAuto.

Calculating the solution¶

And finally, you call:

template<typename Solver, typename System, Number T, template<typename, typename> typename S>
index_t solve_ivp(const ODEProblem<System>& ode, Vector<T, BackendLocation, S>& x, Solver& solver);

to solve the ODEProblem you created in the previous step. Upon entry, x contains the initial conditions of the problem. The value of x is then overwritten with the solution of the problem at the end of the procedure.

Note

You can reuse the same solver for different calls to solve_ivp()

Currently, HighFM supports the following solvers:

template<Number T, DataLocation L = Host> class ForwardEuler¶: The (forward) Euler scheme is the simplest method for solving a system of differential equations. It has a local order of \( O(h^2)\) and a global order of \( O(h)\), where \(h\) is the time step. The value of \(h\) is maintained constant throughout the entire algorithm. Note that the algorithm may diverge if the value of h is too high.

template<Number T, DataLocation L = Host> class RungeKutta4¶: The classical Runge-Kutta of 4th order ( \( O(h^4)\), where \(h\) is the time step). The value of \(h\) is maintained constant throughout the entire algorithm. Note that the algorithm may diverge if the value of h is too high.

template<Number T, DataLocation L = Host> using HighFM::DormandPrince54 = AdaptiveRungeKutta<T, L, internal::DP54Tableau<T>>¶

The adaptive Runge-Kutta method proposed by Dormand and Prince [1]. It calculates the solution with a precision of \( O(h^5)\) and estimates the error of the current step using an \( O(h^4)\) scheme, where \(h\) is the time step. The algorithm then compares the error with a given tolerance and adjusts the time step accordingly.

References:

[1] J. R. Dormand and P. J. Prince, ‘A family of embedded Runge-Kutta formulae’, Journal of Computational and Applied Mathematics, vol. 6, no. 1, pp. 19–26, Mar. 1980, doi: 10.1016/0771-050X(80)90013-3.

template<Number T, DataLocation L = Host> using HighFM::CashKarp54 = AdaptiveRungeKutta<T, L, internal::CK54Tableau<T>>¶

The adaptive Runge-Kutta method proposed by Cash and Karp [1]. It calculates the solution with a precision of \( O(h^5)\) and estimates the error of the current step using an \( O(h^4)\) scheme, where \(h\) is the time step. The algorithm then compares the error with a given tolerance and adjusts the time step accordingly.

References:

[1] J. R. Cash and A. H. Karp, “A variable order Runge-Kutta method for initial value problems with

rapidly varying right-hand sides”, ACM Trans. Math. Softw., vol. 16, no. 3, pp. 201–222, Sep. 1990, doi: 10.1145/79505.79507.