ODE Solvers¶
HighFM support several methods for solving systems of differential equations that have the following form
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 builtin 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) \]

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) \]\[\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 ()
with the following interface:
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 righthand 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 onatol
andrtol
.

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 userdefined system (such that calling
system()
will evalute the righthand side for a given time t).

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

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

std::array<double, 2> time_span¶
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 RungeKutta 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 RungeKutta method proposed by DormandPrince [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 RungeKutta formulae’, Journal of Computational and Applied Mathematics, vol. 6, no. 1, pp. 19–26, Mar. 1980, doi: 10.1016/0771050X(80)900133.

template<Number T, DataLocation L = Host>
using HighFM::CashKarp54 = AdaptiveRungeKutta<T, L, internal::CK54Tableau<T>>¶ The adaptive RungeKutta method proposed by CashKarp [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 RungeKutta method for initial value problems with
rapidly varying righthand sides”, ACM Trans. Math. Softw., vol. 16, no. 3, pp. 201–222, Sep. 1990, doi: 10.1145/79505.79507.