Sandals  v0.0.0
A C++ library for ODEs/DAEs integration
Loading...
Searching...
No Matches
Public Types | Public Member Functions | Private Types | Private Attributes | List of all members
Sandals::RungeKutta< S, N, M > Class Template Reference

Class container for the generic implicit, explicit, and diagonally implicit Runge-Kutta methods. More...

#include <RungeKutta.hxx>

Public Types

using System = typename Implicit<N, M>::Pointer
 
using Type = typename Tableau<S>::Type
 
using Time = Eigen::Vector<Real, Eigen::Dynamic>
 

Public Member Functions

 RungeKutta (const RungeKutta &)=delete
 
RungeKutta & operator= (RungeKutta const &)=delete
 
 RungeKutta (Tableau< S > const &t_tableau)
 
 RungeKutta (Tableau< S > const &t_tableau, System t_system)
 
Type type () const
 
bool is_erk () const
 
bool is_irk () const
 
bool is_dirk () const
 
Tableau< S > & tableau ()
 
Tableau< S > const & tableau () const
 
Integer stages () const
 
std::string name () const
 
Integer order () const
 
bool is_embedded () const
 
MatrixS A () const
 
VectorS b () const
 
VectorS b_embedded () const
 
VectorS c () const
 
System system ()
 
void system (System t_system)
 
bool has_system ()
 
Real absolute_tolerance ()
 
void absolute_tolerance (Real t_absolute_tolerance)
 
Real relative_tolerance ()
 
void relative_tolerance (Real t_relative_tolerance)
 
Real & safety_factor ()
 
void safety_factor (Real t_safety_factor)
 
Real & min_safety_factor ()
 
void min_safety_factor (Real t_min_safety_factor)
 
Real & max_safety_factor ()
 
void max_safety_factor (Real t_max_safety_factor)
 
Real & min_step ()
 
void min_step (Real t_min_step)
 
Integer & max_substeps ()
 
void max_substeps (Integer t_max_substeps)
 
bool adaptive_mode ()
 
void adaptive (bool t_adaptive)
 
void enable_adaptive_mode ()
 
void disable_adaptive_mode ()
 
bool verbose_mode ()
 
void verbose_mode (bool t_verbose)
 
void enable_verbose_mode ()
 
void disable_verbose_mode ()
 
bool reverse_mode ()
 
void reverse (bool t_reverse)
 
void enable_reverse_mode ()
 
void disable_reverse_mode ()
 
Real projection_tolerance ()
 
void projection_tolerance (Real t_projection_tolerance)
 
Integer & max_projection_iterations ()
 
void max_projection_iterations (Integer t_max_projection_iterations)
 
bool projection ()
 
void projection (bool t_projection)
 
void enable_projection ()
 
void disable_projection ()
 
Real estimate_step (VectorN const &x, VectorN const &x_e, Real h_k) const
 
std::string info () const
 
void info (std::ostream &os)
 
bool erk_explicit_step (VectorN const &x_old, Real t_old, Real h_old, VectorN &x_new, Real &h_new) const
 
void erk_implicit_function (Integer s, VectorN const &x, Real t, Real h, MatrixK const &K, VectorN &fun) const
 
void erk_implicit_jacobian (Integer s, VectorN const &x, Real t, Real h, MatrixK const &K, MatrixN &jac) const
 
bool erk_implicit_step (VectorN const &x_old, Real t_old, Real h_old, VectorN &x_new, Real &h_new)
 
void irk_function (VectorN const &x, Real t, Real h, VectorK const &K, VectorK &fun) const
 
void irk_jacobian (VectorN const &x, Real t, Real h, VectorK const &K, MatrixJ &jac) const
 
bool irk_step (VectorN const &x_old, Real t_old, Real h_old, VectorN &x_new, Real &h_new)
 
void dirk_function (Integer n, VectorN const &x, Real t, Real h, MatrixK const &K, VectorN &fun) const
 
void dirk_jacobian (Integer n, VectorN const &x, Real t, Real h, MatrixK const &K, MatrixN &jac) const
 
bool dirk_step (VectorN const &x_old, Real t_old, Real h_old, VectorN &x_new, Real &h_new)
 
bool step (VectorN const &x_old, Real t_old, Real h_old, VectorN &x_new, Real &h_new)
 
bool advance (VectorN const &x_old, Real t_old, Real h_old, VectorN &x_new, Real &h_new)
 
bool solve (VectorX const &t_mesh, VectorN const &ics, Solution< N, M > &sol)
 
bool adaptive_solve (VectorX const &t_mesh, VectorN const &ics, Solution< N, M > &sol)
 
bool project (VectorN const &x, Real t, VectorN &x_projected)
 
bool project_ics (VectorN const &x, Real t, std::vector< Integer > const &projected_equations, std::vector< Integer > const &projected_invariants, VectorN &x_projected) const
 
Real estimate_order (std::vector< VectorX > const &t_mesh, VectorN const &ics, std::function< MatrixX(VectorX)> &sol)
 

Private Types

using VectorK = Eigen::Vector<Real, N*S>
 
using MatrixK = Eigen::Matrix<Real, N, S>
 
using MatrixJ = Eigen::Matrix<Real, N*S, N*S>
 
using VectorP = Eigen::Matrix<Real, N+M, 1>
 
using MatrixP = Eigen::Matrix<Real, N+M, N+M>
 
using NewtonX = Optimist::RootFinder::Newton<N>
 
using NewtonK = Optimist::RootFinder::Newton<N*S>
 
using VectorS = typename Tableau<S>::Vector
 
using MatrixS = typename Tableau<S>::Matrix
 
using VectorN = typename Implicit<N, M>::VectorF
 
using MatrixN = typename Implicit<N, M>::MatrixJF
 
using VectorM = typename Implicit<N, M>::VectorH
 
using MatrixM = typename Implicit<N, M>::MatrixJH
 

Private Attributes

NewtonX m_newtonX
 
NewtonK m_newtonK
 
Tableau< S > m_tableau
 
System m_system
 
Real m_absolute_tolerance {EPSILON_HIGH}
 
Real m_relative_tolerance {EPSILON_HIGH}
 
Real m_safety_factor {0.9}
 
Real m_min_safety_factor {0.1}
 
Real m_max_safety_factor {10.0}
 
Real m_min_step {EPSILON_HIGH}
 
Integer m_max_substeps {5}
 
bool m_adaptive {true}
 
bool m_verbose {false}
 
bool m_reverse {false}
 
Eigen::FullPivLU< MatrixP > m_lu
 
Real m_projection_tolerance {EPSILON_HIGH}
 
Integer m_max_projection_iterations {5}
 
bool m_projection {true}
 

Detailed Description

template<Integer S, Integer N, Integer M>
class Sandals::RungeKutta< S, N, M >

Runge-Kutta methods and Butcher tableaus

In this section, we introduce Runge-Kutta methods for solving ordinary differential equations (ODEs) and differential-algebraic equations (DAEs). Such equation systems describe the dynamics of systems and are used to model a wide range of physical phenomena.

Dynamic systems representation

ODEs and DAEs can be expressed in several forms. While these forms may be mathematically equivalent and convertible into one another, their numerical properties differ, making some more practical for specific computations. Below are the most common representations:

  • Explicit form: \(\mathbf{x}^\prime = \mathbf{f}(\mathbf{x}, t)\).
  • Implicit form: \(\mathbf{F}(\mathbf{x}, \mathbf{x}^\prime, t) = \mathbf{0}\).
  • SemiExplicit form: \(\mathbf{A}(\mathbf{x}, t)\mathbf{x}^\prime = \mathbf{b}(\mathbf{x}, t)\).
  • Linear form: \(\mathbf{E}(t)\mathbf{x}^\prime = \mathbf{A}(t)\mathbf{x} + \mathbf{b}(t)\).

The implicit form is the most general and is frequently used to represent DAEs. However, the explicit forms is more commonly employed for ODEs. However, the explicit, semi-explicit, and linear forms are more commonly employed for ODEs. Note that the semi-explicit and the linear forms are a subset of the explicit form. Consequently, the Runge-Kutta methods we implement will be specialized only for the explicit and implicit forms.

Explicit Runge-Kutta methods

Explicit Runge-Kutta (ERK) methods are among the simplest numerical techniques for solving ordinary differential equations. These methods are termed "explicit" because they typically avoid solving nonlinear systems at each step. Instead, the solution is computed directly using function evaluations. However, this holds only if the derivatives of the dynamic system's states, \(\mathbf{x}^\prime\), can be computed directly or via matrix inversion. This is possible when the dynamic system can be expressed in an explicit, semi-explicit, or linear form.

ERK methods for explicit dynamic systems

For an explicit system of the form \(\mathbf{x}^\prime = \mathbf{f}(\mathbf{x}, t)\), the ERK method is expressed as

\[ \begin{array}{l} \mathbf{K}_i = h_k \mathbf{f} \left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{ij} \mathbf{K}_j, t_k + h_k c_i\right) \\ \mathbf{x}_{k+1} = \mathbf{x}_k + h_k \displaystyle\sum_{i=1}^s b_i \mathbf{K}_i \end{array} \text{.} \]

where \(s\) is the number of stages, \(h_k\) is the step size, \(\mathbf{x}_k\) is the state at time \(t_k\), and \(\mathbf{K}_i\) are the intermediate variables.

For explicit methods, the Runge-Kutta coefficient matrix \(\mathbf{A}\) is strictly lower triangular. Thus, the stages can be computed sequentially as

\[ \mathbf{K}_i = h_k \mathbf{f} \left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{i-1} a_{ij} \mathbf{K}_j, t_k + h_k c_i\right) \text{,} \]

and by unrolling the stages, we obtain

\[ \begin{array}{l} \begin{cases} \mathbf{K}_1 = h_k mathbf{f} \left(\mathbf{x}_k, t_k\right) \\ \mathbf{K}_2 = h_k \mathbf{f} \left(\mathbf{x}_k + a_{21} \mathbf{K}_1, t_k + h_k c_2\right) \\ \vdots \\[-0.5em] \mathbf{K}_s = h_k \mathbf{f} \left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s-1} a_{sj} \mathbf{K}_j, t_k + h_k c_s\right) \end{cases} \\ \mathbf{x}_{k+1} = \mathbf{x}_k + h_k \displaystyle\sum_{i=1}^s b_i \mathbf{K}_i \end{array} \text{.} \]

ERK methods for implicit dynamic systems

For an implicit dynamical system of the form \(\mathbf{F}(\mathbf{x}, \mathbf{x}^\prime, t) = \mathbf{0}\), the ERK method becomes

\[ \begin{array}{l} \mathbf{F}_i \left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{i-1} a_{ij} \mathbf{K}_j, \displaystyle\frac{\mathbf{K}_i}{h_k}, t_k + h_k c_i\right) = \mathbf{0} \\ \mathbf{x}_{k+1} = \mathbf{x}_k + \displaystyle\sum_{i=1}^s b_i \mathbf{K}_i \end{array} \text{.} \]

Here, the \(s\)-stage system of equations \(\mathbf{F}_i\) forms a lower triangular system, typically nonlinear in the intermediate variables \(\mathbf{K}_i\). The solution is obtained using forward substitution, where each stage is solved sequentially, i.e.,

\[ \begin{array}{l} \begin{cases} \mathbf{F}_1 \left(\mathbf{x}_k, \mathbf{K}_1, t_k + h_k c_1\right) = \mathbf{0} \\ \mathbf{F}_2 \left(\mathbf{x}_k + a_{21} \mathbf{K}_1, \displaystyle\frac{\mathbf{K}_2}{h_k}, t_k + h_k c_2\right) = \mathbf{0} \\ \vdots \\[-0.5em] \mathbf{F}_s \left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s-1} a_{sj} \mathbf{K}_j, \displaystyle\frac{\mathbf{K}_s}{h_k}, t_k + h_k c_s\right) = \mathbf{0} \end{cases} \\ \mathbf{x}_{k+1} = \mathbf{x}_k + h_k \displaystyle\sum_{i=1}^s b_i \mathbf{K}_i \end{array} \text{.} \]

Notice that the intermediate variables \(\mathbf{K}_i\) are typically computed using Newton's method.

Implicit Runge-Kutta methods

Implicit Runge-Kutta (IRK) methods are more general than ERK methods, as they can handle implicit systems of the form \(\mathbf{F}(\mathbf{x}, \mathbf{x}^\prime, t) = \mathbf{0}\). IRK methods are more stable but computationally expensive because they require solving a nonlinear system at each step. A key characteristic of IRK methods is that the Runge-Kutta coefficient matrix \(\mathbf{A}\) is fully populated, meaning that the stages are coupled and must be solved simultaneously.

IRK methods for explicit dynamic systems

For an explicit dynamical system of the form \(\mathbf{x}^\prime = \mathbf{f}(\mathbf{x}, t)\), the IRK method is expressed as

\[ \begin{array}{l} \mathbf{K}_i = h_k \mathbf{f} \left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{ij} \mathbf{K}_j, t_k + h_k c_i\right) \\ \mathbf{x}_{k+1} = \mathbf{x}_k + h_k \displaystyle\sum_{i=1}^s b_i \mathbf{K}_i \end{array} \text{,} \]

Unrolling the stages gives

\[ \begin{array}{l} \begin{cases} \mathbf{K}_1 = h_k \mathbf{f} \left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{1j} \mathbf{K}_j, t_k + h_k c_1\right) \\ \mathbf{K}_2 = h_k \mathbf{f} \left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{2j} \mathbf{K}_j, t_k + h_k c_2\right) \\[-0.5em] \vdots \\[-0.5em] \mathbf{K}_s = h_k \mathbf{f} \left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{sj} \mathbf{K}_j, t_k + h_k c_s\right) \end{cases} \\ \mathbf{x}_{k+1} = \mathbf{x}_k + h_k \displaystyle\sum_{i=1}^s b_i \mathbf{K}_i \end{array} \text{.} \]

IRK methods for implicit dynamic systems

For an implicit dynamic system of the form \(\mathbf{F}(\mathbf{x}, \mathbf{x}^\prime, t) = \mathbf{0}\), the IRK method is expressed as

\[ \begin{array}{l} \mathbf{F}_i \left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{ij} \mathbf{K}_j, \displaystyle\frac{\mathbf{K}_i}{h_k}, t_k + h_k c_i\right) = \mathbf{0} \\ \mathbf{x}_{k+1} = \mathbf{x}_k + h_k \displaystyle\sum_{i=1}^s b_i \mathbf{K}_i \end{array} \text{.} \]

Unrolling the stages, we have

\[ \begin{array}{l} \begin{cases} \mathbf{F}_1 \left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{1j} \mathbf{K}_j, \displaystyle\frac{\mathbf{K}_1}{h_k}, t_k + h_k c_1\right) = \mathbf{0} \\ \mathbf{F}_2 \left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{2j} \mathbf{K}_j, \displaystyle\frac{\mathbf{K}_2}{h_k}, t_k + h_k c_2\right) = \mathbf{0} \\[-0.5em] \vdots \\[-0.5em] \mathbf{F}_s \left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{sj} \mathbf{K}_j, \displaystyle\frac{\mathbf{K}_s}{h_k}, t_k + h_k c_s\right) = \mathbf{0} \end{cases} \\ \mathbf{x}_{k+1} = \mathbf{x}_k + h_k \displaystyle\sum_{i=1}^s b_i \mathbf{K}_i \end{array} \text{.} \]

NOTE: In Sandals, the implementation of IRK methods is unified for both explicit and implicit dynamical systems. IRK methods inherently require solving nonlinear systems, irrespective of the system's form. Creating a specialized version for explicit systems would not yield significant computational advantages, as the complexity of solving the nonlinear system remains. Consequently, only the IRK method for implicit dynamical systems is implemented.

Diagonally implicit Runge-Kutta methods

Diagonally implicit Runge-Kutta (DIRK) methods are a specialized subset of IRK methods. They are characterized by a lower triangular \(\mathbf{A}\)-matrix with at least one nonzero diagonal entry. This structure allows each stage to be solved sequentially rather than simultaneously, making DIRK methods more computationally efficient compared to general IRK methods.

DIRK methods for explicit dynamic systems

For an explicit dynamic system of the form \(\mathbf{x}^\prime = \mathbf{f}(\mathbf{x}, t)\), the DIRK method is expressed as

\[ \begin{array}{l} \mathbf{K}_i = h_k \mathbf{f} \left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{i} a_{ij} \mathbf{K}_j, t_k + h_k c_i\right) \\ \mathbf{x}_{k+1} = \mathbf{x}_k + h_k \displaystyle\sum_{i=1}^s b_i \mathbf{K}_i \end{array} \text{.} \]

When the stages are unrolled, the equations become

\[ \begin{array}{l} \begin{cases} \mathbf{K}_1 = h_k \mathbf{f} \left(\mathbf{x}_k + a_{11} \mathbf{K}_1, t_k + h_k c_1\right) \\ \mathbf{K}_2 = h_k \mathbf{f} \left(\mathbf{x}_k + a_{21} \mathbf{K}_1 + h_k a_{22} \mathbf{K}_2, t_k + h_k c_2\right) \\ \vdots \\[-0.5em] \mathbf{K}_s = h_k \mathbf{f} \left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{sj} \mathbf{K}_j, t_k + h_k c_s\right) \end{cases} \\ \mathbf{x}_{k+1} = \mathbf{x}_k + h_k \displaystyle\sum_{i=1}^s b_i \mathbf{K}_i \end{array} \text{.} \]

DIRK methods for implicit dynamic systems

For an implicit dynamic system of the form \(\mathbf{F}(\mathbf{x}, \mathbf{x}^\prime, t) = \mathbf{0}\), the DIRK method is written as

\[ \begin{array}{l} \mathbf{F}_i \left( \mathbf{x}_k + \displaystyle\sum_{j=1}^{i} a_{ij} \mathbf{K}_j, \displaystyle\frac{\mathbf{K}_i}{h_k}, t_k + h_k c_i \right) = \mathbf{0} \\ \mathbf{x}_{k+1} = \mathbf{x}_k + h_k \displaystyle\sum_{i=1}^s b_i \mathbf{K}_i \end{array} \text{.} \]

If the stages are unrolled, we have

\[ \begin{array}{l} \begin{cases} \mathbf{F}_1 \left(\mathbf{x}_k + a_{11} \mathbf{K}_1, \displaystyle\frac{\mathbf{K}_1}{h_k}, t_k + h_k c_1\right) = \mathbf{0} \\ \mathbf{F}_2 \left(\mathbf{x}_k + a_{21} \mathbf{K}_1 + a_{22} \mathbf{K}_2, \displaystyle\frac{\mathbf{K}_3}{h_k}, t_k + h_k c_2\right) = \mathbf{0} \\ \vdots \\[-0.5em] \mathbf{F}_s \left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{sj} \mathbf{K}_j, \displaystyle\frac{\mathbf{K}_s}{h_k}, t_k + h_k c_s\right) = \mathbf{0} \end{cases} \\ \mathbf{x}_{k+1} = \mathbf{x}_k + h_k \displaystyle\sum_{i=1}^s b_i \mathbf{K}_i \end{array} \text{.} \]

NOTE: As with IRK methods, the implementation of DIRK methods for explicit and implicit dynamic systems in Sandals is unified. Since solving the nonlinear system is inherent to DIRK methods, a specialized version for explicit systems would not significantly reduce computational complexity. Consequently, only the DIRK method for implicit dynamic systems is implemented.

Projection on the invariants manifold

In many applications, the dynamic system's states must satisfy certain constraints, which are typically expressed as invariants. For example, in mechanical systems, the constraints may represent the conservation of energy or momentum. In such cases, the states must be projected onto the manifold defined by the invariants to ensure that the system's dynamics are physically consistent. In Sandals, the projection is performed using a simple iterative method, where the states are updated until the invariants are satisfied within a specified tolerance. In other words, at each step, the states are adjusted to ensure that the invariants are preserved.

Constrained minimization problem

Consider the constrained minimization problem

\[ \underset{\mathbf{x}}{\text{minimize}} \quad \frac{1}{2}\left(\mathbf{x} - \tilde{\mathbf{x}}\right)^2 \quad \text{subject to} \quad \mathbf{h}(\mathbf{x}) = \mathbf{0} \text{,} \]

whose Lagrangian is given by

\[ \mathcal{L}(\mathbf{x}, \boldsymbol{\lambda}) = \frac{1}{2}\left(\mathbf{x} - \tilde{\mathbf{x}}\right)^2 + \boldsymbol{\lambda} \cdot \mathbf{h}(\mathbf{x}) \text{.} \]

The first-order Karush-Kuhn-Tucker conditions for this problem are then

\[ \begin{cases} \mathbf{x} + \mathbf{Jh}_{\mathbf{x}}^\top \boldsymbol{\lambda} = \tilde{\mathbf{x}} \\ \mathbf{h}(\mathbf{x}) = \mathbf{0} \end{cases} \text{.} \]

This system of equations can be solved though the iterative solution of a linear system derived from the Taylor expansion

\[ \begin{cases} \mathbf{x} + \delta\mathbf{x} + \mathbf{Jh}_{\mathbf{x}}^\top(\mathbf{x} + \delta\mathbf{x}, t) \boldsymbol{\lambda} = \tilde{\mathbf{x}} \\ \mathbf{h}(\mathbf{x}) + \mathbf{Jh}_{\mathbf{x}}(\mathbf{x}, t) \delta\mathbf{x} + \mathcal{O}\left(\| \delta\mathbf{x} \|^2\right) = \mathbf{0} \end{cases} \text{,} \]

where \(\mathbf{Jh}_{\mathbf{x}}\) is the Jacobian of the constraint function \(\mathbf{h}(\mathbf{x})\) with respect to the states \(\mathbf{x}\). The linear system can be written in matrix form as

\[ \begin{bmatrix} \mathbf{I} & \mathbf{Jh}_{\mathbf{x}}^\top \\ \mathbf{Jh}_{\mathbf{x}} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \delta\mathbf{x} \\ \boldsymbol{\lambda} \end{bmatrix} = \begin{bmatrix} \tilde{\mathbf{x}} - \mathbf{x} \\ -\mathbf{h}(\mathbf{x}) \end{bmatrix} \text{,} \]

and the update step for the states is then \(\mathbf{x} = \tilde{\mathbf{x}} + \delta\mathbf{x}\).

Template Parameters
SThe number of stages of the Runge-Kutta method.
NThe dimension of the ODE/DAE system.
MThe dimension of the invariants manifold.

Member Typedef Documentation

◆ MatrixJ

template<Integer S, Integer N, Integer M>
using Sandals::RungeKutta< S, N, M >::MatrixJ = Eigen::Matrix<Real, N*S, N*S>
private

Templetized matrix type.

◆ MatrixK

template<Integer S, Integer N, Integer M>
using Sandals::RungeKutta< S, N, M >::MatrixK = Eigen::Matrix<Real, N, S>
private

Templetized matrix type.

◆ MatrixM

template<Integer S, Integer N, Integer M>
using Sandals::RungeKutta< S, N, M >::MatrixM = typename Implicit<N, M>::MatrixJH
private

Templetized matrix type.

◆ MatrixN

template<Integer S, Integer N, Integer M>
using Sandals::RungeKutta< S, N, M >::MatrixN = typename Implicit<N, M>::MatrixJF
private

Templetized matrix type.

◆ MatrixP

template<Integer S, Integer N, Integer M>
using Sandals::RungeKutta< S, N, M >::MatrixP = Eigen::Matrix<Real, N+M, N+M>
private

Templetized matrix type.

◆ MatrixS

template<Integer S, Integer N, Integer M>
using Sandals::RungeKutta< S, N, M >::MatrixS = typename Tableau<S>::Matrix
private

Templetized matrix type.

◆ NewtonK

template<Integer S, Integer N, Integer M>
using Sandals::RungeKutta< S, N, M >::NewtonK = Optimist::RootFinder::Newton<N*S>
private

Templetized Newton solver for IRK methods.

◆ NewtonX

template<Integer S, Integer N, Integer M>
using Sandals::RungeKutta< S, N, M >::NewtonX = Optimist::RootFinder::Newton<N>
private

Templetized Newton solver for ERK and DIRK methods.

◆ System

template<Integer S, Integer N, Integer M>
using Sandals::RungeKutta< S, N, M >::System = typename Implicit<N, M>::Pointer

Shared pointer to an implicit ODE/DAE system.

◆ Time

template<Integer S, Integer N, Integer M>
using Sandals::RungeKutta< S, N, M >::Time = Eigen::Vector<Real, Eigen::Dynamic>

Templetized vector type for the independent variable (or time).

◆ Type

template<Integer S, Integer N, Integer M>
using Sandals::RungeKutta< S, N, M >::Type = typename Tableau<S>::Type

Runge-Kutta type enumeration.

◆ VectorK

template<Integer S, Integer N, Integer M>
using Sandals::RungeKutta< S, N, M >::VectorK = Eigen::Vector<Real, N*S>
private

Templetized vector type.

◆ VectorM

template<Integer S, Integer N, Integer M>
using Sandals::RungeKutta< S, N, M >::VectorM = typename Implicit<N, M>::VectorH
private

Templetized vector type.

◆ VectorN

template<Integer S, Integer N, Integer M>
using Sandals::RungeKutta< S, N, M >::VectorN = typename Implicit<N, M>::VectorF
private

Templetized vector type.

◆ VectorP

template<Integer S, Integer N, Integer M>
using Sandals::RungeKutta< S, N, M >::VectorP = Eigen::Matrix<Real, N+M, 1>
private

Templetized vector type.

◆ VectorS

template<Integer S, Integer N, Integer M>
using Sandals::RungeKutta< S, N, M >::VectorS = typename Tableau<S>::Vector
private

Templetized vector type.

Constructor & Destructor Documentation

◆ RungeKutta() [1/3]

template<Integer S, Integer N, Integer M>
Sandals::RungeKutta< S, N, M >::RungeKutta ( const RungeKutta< S, N, M > & )
delete

Copy constructor for the timer.

◆ RungeKutta() [2/3]

template<Integer S, Integer N, Integer M>
Sandals::RungeKutta< S, N, M >::RungeKutta ( Tableau< S > const & t_tableau)
inline

Class constructor for the Runge-Kutta method.

Parameters
[in]t_tableauThe Tableau reference.

◆ RungeKutta() [3/3]

template<Integer S, Integer N, Integer M>
Sandals::RungeKutta< S, N, M >::RungeKutta ( Tableau< S > const & t_tableau,
System t_system )
inline

Class constructor for the Runge-Kutta method.

Parameters
[in]t_tableauThe Tableau reference.
[in]t_systemThe ODE/DAE system shared pointer.

Member Function Documentation

◆ A()

template<Integer S, Integer N, Integer M>
MatrixS Sandals::RungeKutta< S, N, M >::A ( ) const
inline

Get the Butcher tableau matrix \( \mathbf{A} \).

Returns
The Butcher tableau matrix \( \mathbf{A} \).

◆ absolute_tolerance() [1/2]

template<Integer S, Integer N, Integer M>
Real Sandals::RungeKutta< S, N, M >::absolute_tolerance ( )
inline

Get the adaptive step absolute tolerance \( \epsilon_{\text{abs}} \).

Returns
The adaptive step absolute tolerance \( \epsilon_{\text{abs}} \).

◆ absolute_tolerance() [2/2]

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::absolute_tolerance ( Real t_absolute_tolerance)
inline

Get the adaptive step absolute tolerance reference \( \epsilon_{\text{abs}} \).

Parameters
[in]t_absolute_toleranceThe adaptive step absolute tolerance \( \epsilon_{\text{abs}} \).

◆ adaptive()

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::adaptive ( bool t_adaptive)
inline

Set the adaptive step mode.

Parameters
[in]t_adaptiveThe adaptive step mode.

◆ adaptive_mode()

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::adaptive_mode ( )
inline

Get the adaptive step mode.

Returns
The adaptive step mode.

◆ adaptive_solve()

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::adaptive_solve ( VectorX const & t_mesh,
VectorN const & ics,
Solution< N, M > & sol )
inline

Solve the system and calculate the approximate solution over the suggested independent variable mesh \( \mathbf{t} = \left[ t_1, t_2, \ldots, t_n \right]^\top \), the step size is automatically computed based on the error control method.

Parameters
[in]t_meshIndependent variable (or time) mesh \( \mathbf{t} \).
[in]icsInitial conditions \( \mathbf{x}(t = 0) \).
[out]solThe solution of the system over the mesh of independent variable.
Returns
True if the system is successfully solved, false otherwise.

◆ advance()

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::advance ( VectorN const & x_old,
Real t_old,
Real h_old,
VectorN & x_new,
Real & h_new )
inline

Advance using a generic integration method for a system of the form \( \mathbf{F}(\mathbf{x}, \mathbf{x}^\prime, t) = \mathbf{0} \). The step is automatically selected based on the system properties and the integration method properties. In the advvancing step, the system solution is also projected on the manifold \( \mathbf{h}(\mathbf{x}, t) \). Substepping may be also used to if the integration step fails with the current step size.

Parameters
[in]x_oldStates \( \mathbf{x}_k \) at the \( k \)-th step.
[in]t_oldIndependent variable (or time) \( t_k \) at the \( k \)-th step.
[in]h_oldAdvancing step \( h_k \) at the \( k \)-th step.
[out]x_newComputed states \( \mathbf{x}_{k+1} \) at the \( (k+1) \)-th step.
[out]h_newThe suggested step \( h_{k+1}^\star \) for the next advancing step.
Returns
True if the step is successfully computed, false otherwise.

◆ b()

template<Integer S, Integer N, Integer M>
VectorS Sandals::RungeKutta< S, N, M >::b ( ) const
inline

Get the Butcher tableau weights vector \( \mathbf{b} \).

Returns
The Butcher tableau weights vector \( \mathbf{b} \).

◆ b_embedded()

template<Integer S, Integer N, Integer M>
VectorS Sandals::RungeKutta< S, N, M >::b_embedded ( ) const
inline

Get the Butcher tableau embedded weights vector \( \hat{\mathbf{b}} \).

Returns
The Butcher tableau embedded weights vector \( \hat{\mathbf{b}} \).

◆ c()

template<Integer S, Integer N, Integer M>
VectorS Sandals::RungeKutta< S, N, M >::c ( ) const
inline

Get the Butcher tableau nodes vector \( \mathbf{c} \).

Returns
The Butcher tableau nodes vector \( \mathbf{c} \).

◆ dirk_function()

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::dirk_function ( Integer n,
VectorN const & x,
Real t,
Real h,
MatrixK const & K,
VectorN & fun ) const
inline

Compute the residual of system to be solved, which is given by the values of the system

\[\begin{array}{l} \mathbf{F}_n \left(\mathbf{x} + \displaystyle\sum_{j=1}^n a_{nj}\tilde{\mathbf{K}}_j, \displaystyle\frac{\tilde{\mathbf{K}}_n}{h}, t + h c_n\right) \end{array} \text{.} \]

where \( \tilde{\mathbf{K}} = h \mathbf{K} \).

Note
If \( h \to 0 \), then the first guess to solve the nonlinear system is given by \(\tilde{\mathbf{K}} = \mathbf{0} \).
Parameters
[in]nEquation index \( n \).
[in]xStates \( \mathbf{x} \).
[in]tIndependent variable (or time) \( t \).
[in]hAdvancing step \( h \).
[in]KVariables \( \tilde{\mathbf{K}} = h \mathbf{K} \) of the system to be solved.
[out]funThe residual of system to be solved.

◆ dirk_jacobian()

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::dirk_jacobian ( Integer n,
VectorN const & x,
Real t,
Real h,
MatrixK const & K,
MatrixN & jac ) const
inline

Compute the Jacobian of the system of equations

\[\begin{array}{l} \mathbf{F}_n \left(\mathbf{x} + \displaystyle\sum_{j=1}^n a_{nj}\tilde{\mathbf{K}}_j, \displaystyle\frac{\tilde{\mathbf{K}}_n}{h}, t + h c_n\right) \end{array} \text{.} \]

with respect to the \( \tilde{\mathbf{K}}_n = h \mathbf{K}_n \) variables. The Jacobian matrix of the \( n \)-th equation with respect to the \( \tilde{\mathbf{K}}_n \) variables is given by

\[ \frac{\partial\mathbf{F}_n}{\partial\tilde{\mathbf{K}}_n}\left(\mathbf{x} + \displaystyle \sum_{j=1}^n a_{nj}\tilde{\mathbf{K}}_j, \displaystyle\frac{\tilde{\mathbf{K}}_n}{h}, t + h c_n \right) = a_{nj} \mathbf{JF}_x + \displaystyle\frac{1}{h}\begin{cases} \mathbf{JF}_{x^{\prime}} & n = j \\ 0 & n \neq j \end{cases} \text{,} \]

for \( j = 1, 2, \ldots, s \).

Note
If \( h \to 0 \), then the first guess to solve the nonlinear system is given by \(\tilde{\mathbf{K}} = \mathbf{0} \).
Parameters
[in]nEquation index \( n \).
[in]xStates \( \mathbf{x} \).
[in]tIndependent variable (or time) \( t \).
[in]hAdvancing step \( h \).
[in]KVariables \( h \mathbf{K} \) of the system to be solved.
[out]jacThe Jacobian of system to be solved.

◆ dirk_step()

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::dirk_step ( VectorN const & x_old,
Real t_old,
Real h_old,
VectorN & x_new,
Real & h_new )
inline

Compute the new states \( \mathbf{x}_{k+1} \) at the next advancing step \( t_{k+1} = t_k + h_k \) using an diagonally implicit Runge-Kutta method, given an implicit system of the form \( \mathbf{F}(\mathbf{x}, \mathbf{x}^\prime, t) = \mathbf{0} \). If the method is embedded, the step for the next advancing step \( h_{k+1}^\star \) is also computed according to the error control method.

Parameters
[in]x_oldStates \( \mathbf{x}_k \) at the \( k \)-th step.
[in]t_oldIndependent variable (or time) \( t_k \) at the \( k \)-th step.
[in]h_oldAdvancing step \( h_k \) at the \( k \)-th step.
[out]x_newComputed states \( \mathbf{x}_{k+1} \) at the \( (k+1) \)-th step.
[out]h_newThe suggested step \( h_{k+1}^\star \) for the next advancing step.
Returns
True if the step is successfully computed, false otherwise.

◆ disable_adaptive_mode()

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::disable_adaptive_mode ( )
inline

Disable the adaptive step mode.

◆ disable_projection()

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::disable_projection ( )
inline

Disable the projection mode.

◆ disable_reverse_mode()

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::disable_reverse_mode ( )
inline

Disable the time reverse mode.

◆ disable_verbose_mode()

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::disable_verbose_mode ( )
inline

Disable the verbose mode.

◆ enable_adaptive_mode()

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::enable_adaptive_mode ( )
inline

Enable the adaptive step mode.

◆ enable_projection()

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::enable_projection ( )
inline

Enable the projection mode.

◆ enable_reverse_mode()

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::enable_reverse_mode ( )
inline

Enable the time reverse mode.

◆ enable_verbose_mode()

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::enable_verbose_mode ( )
inline

Enable the verbose mode.

◆ erk_explicit_step()

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::erk_explicit_step ( VectorN const & x_old,
Real t_old,
Real h_old,
VectorN & x_new,
Real & h_new ) const
inline

Compute the new states \( \mathbf{x}_{k+1} \) at the next advancing step \( t_{k+1} = t_k + h_k \) using an explicit Runge-Kutta method, given an explicit system of the form \(\mathbf{x}^\prime = \mathbf{f}(\mathbf{x}, t) \). If the method is embedded, the step for the next advancing step \( h_{k+1}^\star \) is also computed according to the error control method.

Parameters
[in]x_oldStates \( \mathbf{x}_k \) at the \( k \)-th step.
[in]t_oldIndependent variable (or time) \( t_k \) at the \( k \)-th step.
[in]h_oldAdvancing step \( h_k \) at the \( k \)-th step.
[out]x_newComputed states \( \mathbf{x}_{k+1} \) at the \( (k+1) \)-th step.
[out]h_newThe suggested step \( h_{k+1}^\star \) for the next advancing step.
Returns
True if the step is successfully computed, false otherwise.

◆ erk_implicit_function()

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::erk_implicit_function ( Integer s,
VectorN const & x,
Real t,
Real h,
MatrixK const & K,
VectorN & fun ) const
inline

Compute the residual of system to be solved, which is given by the values of the system

\[\begin{array}{l} \mathbf{F}_n \left(\mathbf{x} + \displaystyle\sum_{j=1}^{n-1} a_{nj}\tilde{\mathbf{K}}_j, \displaystyle\frac{\tilde{\mathbf{K}}_n}{h}, t + h c_n\right) \end{array} \text{.} \]

where \( \tilde{\mathbf{K}} = h \mathbf{K} \).

Note
If \( h \to 0 \), then the first guess to solve the nonlinear system is given by \(\tilde{\mathbf{K}} = \mathbf{0} \).
Parameters
[in]sStage index \( s \).
[in]xStates \( \mathbf{x} \).
[in]tIndependent variable (or time) \( t \).
[in]hAdvancing step \( h \).
[in]KVariables \( \tilde{\mathbf{K}} = h \mathbf{K} \) of the system to be solved.
[out]funThe residual of system to be solved.

◆ erk_implicit_jacobian()

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::erk_implicit_jacobian ( Integer s,
VectorN const & x,
Real t,
Real h,
MatrixK const & K,
MatrixN & jac ) const
inline

Compute the Jacobian of the system of equations

\[\begin{array}{l} \mathbf{F}_n \left(\mathbf{x} + \displaystyle\sum_{j=1}^{n-1} a_{nj}\tilde{\mathbf{K}}_j, \displaystyle\frac{\tilde{\mathbf{K}}_n}{h}, t + h c_n\right) \end{array} \text{.} \]

with respect to the \( \tilde{\mathbf{K}}_n = h \mathbf{K}_n \) variables. The Jacobian matrix of the \( n \)-th equation with respect to the \( \tilde{\mathbf{K}}_n \) variables is given by

\[ \frac{\partial\mathbf{F}_n}{\partial\tilde{\mathbf{K}}_n}\left(\mathbf{x} + \displaystyle \sum_{j=1}^{n-1} a_{nj}\tilde{\mathbf{K}}_j, \displaystyle\frac{\tilde{\mathbf{K}}_n}{h}, t + h c_n \right) = a_{nj} \mathbf{JF}_x + \displaystyle\frac{1}{h}\begin{cases} \mathbf{JF}_{x^{\prime}} & n = j \\ 0 & n \neq j \end{cases} \text{,} \]

for \( j = 1, 2, \ldots, s \).

Parameters
[in]sStage index \( s \).
[in]xStates \( \mathbf{x} \).
[in]tIndependent variable (or time) \( t \).
[in]hAdvancing step \( h \).
[in]KVariables \( h \mathbf{K} \) of the system to be solved.
[out]jacThe Jacobian of system to be solved.

◆ erk_implicit_step()

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::erk_implicit_step ( VectorN const & x_old,
Real t_old,
Real h_old,
VectorN & x_new,
Real & h_new )
inline

Compute the new states \( \mathbf{x}_{k+1} \) at the next advancing step \( t_{k+1} = t_k + h_k \) using an explicit Runge-Kutta method, given an explicit system of the form \(\mathbf{x}^\prime = \mathbf{f}(\mathbf{x}, t) \). If the method is embedded, the step for the next advancing step \( h_{k+1}^\star \) is also computed according to the error control method.

Parameters
[in]x_oldStates \( \mathbf{x}_k \) at the \( k \)-th step.
[in]t_oldIndependent variable (or time) \( t_k \) at the \( k \)-th step.
[in]h_oldAdvancing step \( h_k \) at the \( k \)-th step.
[out]x_newComputed states \( \mathbf{x}_{k+1} \) at the \( (k+1) \)-th step.
[out]h_newThe suggested step \( h_{k+1}^\star \) for the next advancing step.
Returns
True if the step is successfully computed, false otherwise.

◆ estimate_order()

template<Integer S, Integer N, Integer M>
Real Sandals::RungeKutta< S, N, M >::estimate_order ( std::vector< VectorX > const & t_mesh,
VectorN const & ics,
std::function< MatrixX(VectorX)> & sol )
inline

Estimate the order of the Runge-Kutta method.

Parameters
[in]t_meshThe vector of time meshes with same initial and final time with fixed step.
[in]icsInitial conditions \( \mathbf{x}(t = 0) \).
[in]solThe analytical solution function.
Returns
The estimated order of the method.

◆ estimate_step()

template<Integer S, Integer N, Integer M>
Real Sandals::RungeKutta< S, N, M >::estimate_step ( VectorN const & x,
VectorN const & x_e,
Real h_k ) const
inline

Estimate the optimal step size for the next advancing step according to the error control method. The error control method used is based on the local truncation error, which is computed as

\[\varepsilon_{\text{t}} = \max\left(\left| \mathbf{x} - \hat{\mathbf{x}} \right|\right) \]

where \( \mathbf{x} \) is the approximation of the states computed with the higher order method, and \( \hat{\mathbf{x}} \) is the approximation of the states computed with the lower order method. The desired error is computed as

\[\varepsilon_{\text{d}} = \epsilon_{\text{abs}} + \epsilon_{\text{rel}} \max\left(\left| \mathbf{x} \right|, \left| \hat{\mathbf{x}} \right|\right) \]

where \( \epsilon_{\text{abs}} \) is the absolute tolerance and \( \epsilon_{\text{rel}} \) is the relative tolerance. The optimal step size is then estimated as

\[h_{\text{opt}} = h_k \min\left( f_{\max}, \max\left( f_{\min}, f \left( \displaystyle\frac{\varepsilon_{\text{d}}}{\varepsilon_{\text{t}}} \right)^{\frac{1}{q+1}} \right)\right) \]

where \( f \) is the safety factor, \( f_{\max} \) is the maximum safety factor, \( f_{\min} \) is the minimum safety factor, and \( q \) is the order of the embedded method.

Parameters
[in]xStates approximation \( \mathbf{x}_{k+1} \).
[in]x_eEmbedded method's states approximation \( \hat{\mathbf{x}}_{k+1} \).
[in]h_kActual advancing step \( h_k \).
Returns
The suggested step for the next integration step \( h_{k+1}^\star \).

◆ has_system()

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::has_system ( )
inline

Check if the ODE/DAE system pointer is set.

Returns
True if the ODE/DAE system pointer is set, false otherwise.

◆ info() [1/2]

template<Integer S, Integer N, Integer M>
std::string Sandals::RungeKutta< S, N, M >::info ( ) const
inline

Print the Runge-Kutta method information to the output stream.

Returns
The Runge-Kutta method information as a string.

◆ info() [2/2]

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::info ( std::ostream & os)
inline

Print the Runge-Kutta method information on a stream.

Parameters
[in,out]osOutput stream.

◆ irk_function()

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::irk_function ( VectorN const & x,
Real t,
Real h,
VectorK const & K,
VectorK & fun ) const
inline

Compute the residual of system to be solved, which is given by the values of the system

\[\begin{cases} \mathbf{F}_1 \left(\mathbf{x} + \displaystyle\sum_{j=1}^{s} a_{1j}\tilde{\mathbf{K}}_j, \displaystyle\frac{\tilde{\mathbf{K}}_1}{h}, t + h c_1\right) \\ \mathbf{F}_2 \left(\mathbf{x} + \displaystyle\sum_{j=1}^{s} a_{2j}\tilde{\mathbf{K}}_j, \displaystyle\frac{\tilde{\mathbf{K}}_2}{h}, t + h c_2\right) \\[-0.5em] \vdots \\[-0.5em] \mathbf{F}_s \left(\mathbf{x} + \displaystyle\sum_{j=1}^{s} a_{sj}\tilde{\mathbf{K}}_j, \displaystyle\frac{\tilde{\mathbf{K}}_s}{h}, t + h c_s\right) \end{cases} \]

where \( \tilde{\mathbf{K}} = h \mathbf{K} \).

Note
If \( h \to 0 \), then the first guess to solve the nonlinear system is given by \(\tilde{\mathbf{K}} = \mathbf{0} \).
Parameters
[in]xStates \( \mathbf{x} \).
[in]tIndependent variable (or time) \( t \).
[in]hAdvancing step \( h \).
[in]KVariables \( \tilde{\mathbf{K}} = h \mathbf{K} \) of the system to be solved.
[out]funThe residual of system to be solved.

◆ irk_jacobian()

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::irk_jacobian ( VectorN const & x,
Real t,
Real h,
VectorK const & K,
MatrixJ & jac ) const
inline

Compute the Jacobian of the system of equations

\[\begin{cases} \mathbf{F}_1 \left(\mathbf{x} + \displaystyle\sum_{j=1}^{s} a_{1j}\tilde{\mathbf{K}}_j, \displaystyle\frac{\tilde{\mathbf{K}}_1}{h}, t + h c_1\right) \\ \mathbf{F}_2 \left(\mathbf{x} + \displaystyle\sum_{j=1}^{s} a_{2j}\tilde{\mathbf{K}}_j, \displaystyle\frac{\tilde{\mathbf{K}}_2}{h}, t + h c_2\right) \\[-0.5em] \vdots \\[-0.5em] \mathbf{F}_s \left(\mathbf{x} + \displaystyle\sum_{j=1}^{s} a_{sj}\tilde{\mathbf{K}}_j, \displaystyle\frac{\tilde{\mathbf{K}}_s}{h}, t + h c_s\right) \end{cases} \text{.} \]

with respect to the \( \tilde{\mathbf{K}} = h \mathbf{K} \) variables. The \( ij \)-th Jacobian matrix entry is given by

\[ \frac{\partial\mathbf{F}_i}{\partial\tilde{\mathbf{K}}_j}\left(\mathbf{x} + \displaystyle \sum_{j=1}^s a_{ij}\tilde{\mathbf{K}}_j, \displaystyle\frac{\tilde{\mathbf{K}}_i}{h}, t + h c_i \right) = a_{ij} \mathbf{JF}_x + \displaystyle\frac{1}{h}\begin{cases} \mathbf{JF}_{x^{\prime}} & i = j \\ 0 & i \neq j \end{cases} \text{,} \]

for \( i = 1, 2, \ldots, s \) and \( j = 1, 2, \ldots, s \).

Note
If \( h \to 0 \), then the first guess to solve the nonlinear system is given by \(\tilde{\mathbf{K}} = \mathbf{0} \).
Parameters
[in]xStates \( \mathbf{x} \).
[in]tIndependent variable (or time) \( t \).
[in]hAdvancing step \( h \).
[in]KVariables \( h \mathbf{K} \) of the system to be solved.
[out]jacThe Jacobian of system to be solved.

◆ irk_step()

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::irk_step ( VectorN const & x_old,
Real t_old,
Real h_old,
VectorN & x_new,
Real & h_new )
inline

Compute the new states \( \mathbf{x}_{k+1} \) at the next advancing step \( t_{k+1} = t_k + h_k \) using an implicit Runge-Kutta method, given an implicit system of the form \(\mathbf{F}(\mathbf{x}, \mathbf{x}^\prime, t) = \mathbf{0} \). If the method is embedded, the step for the next advancing step \( h_{k+1}^\star \) is also computed according to the error control method.

Parameters
[in]x_oldStates \( \mathbf{x}_k \) at the \( k \)-th step.
[in]t_oldIndependent variable (or time) \( t_k \) at the \( k \)-th step.
[in]h_oldAdvancing step \( h_k \) at the \( k \)-th step.
[out]x_newComputed states \( \mathbf{x}_{k+1} \) at the \( (k+1) \)-th step.
[out]h_newThe suggested step \( h_{k+1}^\star \) for the next advancing step.
Returns
True if the step is successfully computed, false otherwise.

◆ is_dirk()

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::is_dirk ( ) const
inline

Check if the method is a diagonally implicit Runge-Kutta (DIRK) method.

Returns
True if the method is a diagonally implicit Runge-Kutta method, false otherwise.

◆ is_embedded()

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::is_embedded ( ) const
inline

Check if the Runge-Kutta method is embedded.

Returns
True if the Runge-Kutta method is embedded, false otherwise.

◆ is_erk()

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::is_erk ( ) const
inline

Check if the method is an explicit Runge-Kutta (ERK) method.

Returns
True if the method is an explicit Runge-Kutta method, false otherwise.

◆ is_irk()

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::is_irk ( ) const
inline

Check if the method is an implicit Runge-Kutta (IRK) method.

Returns
True if the method is an implicit Runge-Kutta method, false otherwise.

◆ max_projection_iterations() [1/2]

template<Integer S, Integer N, Integer M>
Integer & Sandals::RungeKutta< S, N, M >::max_projection_iterations ( )
inline

Get the maximum number of projection iterations.

Returns
The maximum number of projection iterations.

◆ max_projection_iterations() [2/2]

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::max_projection_iterations ( Integer t_max_projection_iterations)
inline

Set the maximum number of projection iterations.

Parameters
[in]t_max_projection_iterationsThe maximum number of projection iterations.

◆ max_safety_factor() [1/2]

template<Integer S, Integer N, Integer M>
Real & Sandals::RungeKutta< S, N, M >::max_safety_factor ( )
inline

Get the maximum safety factor for adaptive step \( f_{\max} \).

Returns
The maximum safety factor for adaptive step \( f_{\max} \).

◆ max_safety_factor() [2/2]

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::max_safety_factor ( Real t_max_safety_factor)
inline

Set the maximum safety factor for adaptive step \( f_{\max} \).

Parameters
[in]t_max_safety_factorThe maximum safety factor for adaptive step \( f_{\max} \).

◆ max_substeps() [1/2]

template<Integer S, Integer N, Integer M>
Integer & Sandals::RungeKutta< S, N, M >::max_substeps ( )
inline

Get the maximum number of substeps.

Returns
The maximum number of substeps.

◆ max_substeps() [2/2]

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::max_substeps ( Integer t_max_substeps)
inline

Set the maximum number of substeps.

Parameters
[in]t_max_substepsThe maximum number of substeps.

◆ min_safety_factor() [1/2]

template<Integer S, Integer N, Integer M>
Real & Sandals::RungeKutta< S, N, M >::min_safety_factor ( )
inline

Get the minimum safety factor for adaptive step \( f_{\min} \).

Returns
The minimum safety factor for adaptive step \( f_{\min} \).

◆ min_safety_factor() [2/2]

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::min_safety_factor ( Real t_min_safety_factor)
inline

Set the minimum safety factor for adaptive step \( f_{\min} \).

Parameters
[in]t_min_safety_factorThe minimum safety factor for adaptive step \( f_{\min} \).

◆ min_step() [1/2]

template<Integer S, Integer N, Integer M>
Real & Sandals::RungeKutta< S, N, M >::min_step ( )
inline

Get the minimum step for advancing \( h_{\min} \).

Returns
The minimum step for advancing \( h_{\min} \).

◆ min_step() [2/2]

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::min_step ( Real t_min_step)
inline

Set the minimum step for advancing \( h_{\min} \).

Parameters
[in]t_min_stepThe minimum step for advancing \( h_{\min} \).

◆ name()

template<Integer S, Integer N, Integer M>
std::string Sandals::RungeKutta< S, N, M >::name ( ) const
inline

Get the name of the Runge-Kutta method.

Returns
The name of the Runge-Kutta method.

◆ operator=()

template<Integer S, Integer N, Integer M>
RungeKutta & Sandals::RungeKutta< S, N, M >::operator= ( RungeKutta< S, N, M > const & )
delete

Assignment operator for the timer.

◆ order()

template<Integer S, Integer N, Integer M>
Integer Sandals::RungeKutta< S, N, M >::order ( ) const
inline

Get the order \( p \) of the Runge-Kutta method.

Returns
The order \( p \) of the Runge-Kutta method.

◆ project()

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::project ( VectorN const & x,
Real t,
VectorN & x_projected )
inline

Project the system solution \( \mathbf{x} \) on the invariants \( \mathbf{h} (\mathbf{x}, t) = \mathbf{0} \).

Parameters
[in]xThe initial guess for the states \(\widetilde{\mathbf{x}} \).
[in]tThe independent variable (or time) \( t \) at which the states are evaluated.
[out]x_projectedThe projected states \( \mathbf{x} \) closest to the invariants manifold \( \mathbf{h} (\mathbf{x}, t) = \mathbf{0} \).
Returns
True if the solution is successfully projected, false otherwise.

◆ project_ics()

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::project_ics ( VectorN const & x,
Real t,
std::vector< Integer > const & projected_equations,
std::vector< Integer > const & projected_invariants,
VectorN & x_projected ) const
inline

Project the system solution \( \mathbf{x} \) on the invariants \( \mathbf{h} (\mathbf{x}, t) = \mathbf{0} \).

Parameters
[in]xThe initial guess for the states \(\widetilde{\mathbf{x}} \).
[in]tThe independent variable (or time) \( t \) at which the states are evaluated.
[in]projected_equationsThe indices of the states to be projected.
[in]projected_invariantsThe indices of the invariants to be projected.
[out]x_projectedThe projected states \( \mathbf{x} \) closest to the invariants manifold \( \mathbf{h} (\mathbf{x}, t) = \mathbf{0} \).
Returns
True if the solution is successfully projected, false otherwise.

◆ projection() [1/2]

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::projection ( )
inline

Get projection mode.

Returns
The projection mode.

◆ projection() [2/2]

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::projection ( bool t_projection)
inline

Set the projection mode.

Parameters
[in]t_projectionThe projection mode.

◆ projection_tolerance() [1/2]

template<Integer S, Integer N, Integer M>
Real Sandals::RungeKutta< S, N, M >::projection_tolerance ( )
inline

Get the projection tolerance.

Returns
The projection tolerance.

◆ projection_tolerance() [2/2]

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::projection_tolerance ( Real t_projection_tolerance)
inline

Set the projection tolerance.

Parameters
[in]t_projection_toleranceThe projection tolerance.

◆ relative_tolerance() [1/2]

template<Integer S, Integer N, Integer M>
Real Sandals::RungeKutta< S, N, M >::relative_tolerance ( )
inline

Get the adaptive step relative tolerance \( \epsilon_{\text{rel}} \).

Returns
The adaptive step relative tolerance \( \epsilon_{\text{rel}} \).

◆ relative_tolerance() [2/2]

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::relative_tolerance ( Real t_relative_tolerance)
inline

Get the adaptive step relative tolerance reference \( \epsilon_{\text{rel}} \).

Parameters
[in]t_relative_toleranceThe adaptive step relative tolerance \( \epsilon_{\text{rel}} \).

◆ reverse()

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::reverse ( bool t_reverse)
inline

Set the time reverse mode.

Parameters
[in]t_reverseThe time reverse mode.

◆ reverse_mode()

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::reverse_mode ( )
inline

Get the time reverse mode.

Returns
The time reverse mode.

◆ safety_factor() [1/2]

template<Integer S, Integer N, Integer M>
Real & Sandals::RungeKutta< S, N, M >::safety_factor ( )
inline

Get the safety factor for adaptive step \( f \).

Returns
The safety factor for adaptive step \( f \).

◆ safety_factor() [2/2]

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::safety_factor ( Real t_safety_factor)
inline

Set safety factor for adaptive step \( f \).

Parameters
[in]t_safety_factorThe safety factor for adaptive step \( f \).

◆ solve()

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::solve ( VectorX const & t_mesh,
VectorN const & ics,
Solution< N, M > & sol )
inline

Solve the system and calculate the approximate solution over the independent variable (or time) mesh \( \mathbf{t} = \left[ t_1, t_2, \ldots, t_n \right]^\top \). The step size is fixed and given by \( h = t_{k+1} - t_k \).

Parameters
[in]t_meshIndependent variable (or time) mesh \( \mathbf{t} \).
[in]icsInitial conditions \( \mathbf{x}(t = 0) \).
[out]solThe solution of the system over the mesh of independent variable.
Returns
True if the system is successfully solved, false otherwise.

◆ stages()

template<Integer S, Integer N, Integer M>
Integer Sandals::RungeKutta< S, N, M >::stages ( ) const
inline

Get the stages \( s \) number of the Runge-Kutta method.

Returns
The stages \( s \) number of the Runge-Kutta method.

◆ step()

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::step ( VectorN const & x_old,
Real t_old,
Real h_old,
VectorN & x_new,
Real & h_new )
inline

Compute a step using a generic integration method for a system of the form \( \mathbf{F}( \mathbf{x}, \mathbf{x}^\prime, t) = \mathbf{0} \). The step is automatically selected based on the system properties and the integration method properties.

Parameters
[in]x_oldStates \( \mathbf{x}_k \) at the \( k \)-th step.
[in]t_oldIndependent variable (or time) \( t_k \) at the \( k \)-th step.
[in]h_oldAdvancing step \( h_k \) at the \( k \)-th step.
[out]x_newComputed states \( \mathbf{x}_{k+1} \) at the \( (k+1) \)-th step.
[out]h_newThe suggested step \( h_{k+1}^\star \) for the next advancing step.
Returns
True if the step is successfully computed, false otherwise.

◆ system() [1/2]

template<Integer S, Integer N, Integer M>
System Sandals::RungeKutta< S, N, M >::system ( )
inline

Get the ODE/DAE system pointer.

Returns
The ODE/DAE system pointer.

◆ system() [2/2]

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::system ( System t_system)
inline

Set the ODE/DAE system pointer.

Parameters
[in]t_systemThe ODE/DAE system pointer.

◆ tableau() [1/2]

template<Integer S, Integer N, Integer M>
Tableau< S > & Sandals::RungeKutta< S, N, M >::tableau ( )
inline

Get the Butcher Tableau reference.

Returns
The Tableau reference.

◆ tableau() [2/2]

template<Integer S, Integer N, Integer M>
Tableau< S > const & Sandals::RungeKutta< S, N, M >::tableau ( ) const
inline

Get the Butcher Tableau const reference.

Returns
The Tableau const reference.

◆ type()

template<Integer S, Integer N, Integer M>
Type Sandals::RungeKutta< S, N, M >::type ( ) const
inline

Get the enumeration type of the Runge-Kutta method.

Returns
The enumeration type of the Runge-Kutta method.

◆ verbose_mode() [1/2]

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::verbose_mode ( )
inline

Get the verbose mode.

Returns
The verbose mode.

◆ verbose_mode() [2/2]

template<Integer S, Integer N, Integer M>
void Sandals::RungeKutta< S, N, M >::verbose_mode ( bool t_verbose)
inline

Set the verbose mode.

Parameters
[in]t_verboseThe verbose mode.

Member Data Documentation

◆ m_absolute_tolerance

template<Integer S, Integer N, Integer M>
Real Sandals::RungeKutta< S, N, M >::m_absolute_tolerance {EPSILON_HIGH}
private

Absolute tolerance for adaptive step \( \epsilon_{\text{abs}} \).

◆ m_adaptive

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::m_adaptive {true}
private

Aadaptive step mode boolean.

◆ m_lu

template<Integer S, Integer N, Integer M>
Eigen::FullPivLU<MatrixP> Sandals::RungeKutta< S, N, M >::m_lu
private

LU decomposition for the projection matrix.

◆ m_max_projection_iterations

template<Integer S, Integer N, Integer M>
Integer Sandals::RungeKutta< S, N, M >::m_max_projection_iterations {5}
private

Maximum number of projection steps.

◆ m_max_safety_factor

template<Integer S, Integer N, Integer M>
Real Sandals::RungeKutta< S, N, M >::m_max_safety_factor {10.0}
private

Maximum safety factor for adaptive step \( f_{\min} \).

◆ m_max_substeps

template<Integer S, Integer N, Integer M>
Integer Sandals::RungeKutta< S, N, M >::m_max_substeps {5}
private

Maximum number of substeps.

◆ m_min_safety_factor

template<Integer S, Integer N, Integer M>
Real Sandals::RungeKutta< S, N, M >::m_min_safety_factor {0.1}
private

Minimum safety factor for adaptive step \( f_{\max} \).

◆ m_min_step

template<Integer S, Integer N, Integer M>
Real Sandals::RungeKutta< S, N, M >::m_min_step {EPSILON_HIGH}
private

Minimum step for advancing \( h_{\min} \).

◆ m_newtonK

template<Integer S, Integer N, Integer M>
NewtonK Sandals::RungeKutta< S, N, M >::m_newtonK
private

Newton solver for IRK methods.

◆ m_newtonX

template<Integer S, Integer N, Integer M>
NewtonX Sandals::RungeKutta< S, N, M >::m_newtonX
private

Newton solver for ERK and DIRK methods.

◆ m_projection

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::m_projection {true}
private

Aadaptive step mode boolean.

◆ m_projection_tolerance

template<Integer S, Integer N, Integer M>
Real Sandals::RungeKutta< S, N, M >::m_projection_tolerance {EPSILON_HIGH}
private

Projection tolerance \( \epsilon_{\text{proj}} \).

◆ m_relative_tolerance

template<Integer S, Integer N, Integer M>
Real Sandals::RungeKutta< S, N, M >::m_relative_tolerance {EPSILON_HIGH}
private

Relative tolerance for adaptive step \( \epsilon_{\text{rel}} \).

◆ m_reverse

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::m_reverse {false}
private

Time reverse mode boolean.

◆ m_safety_factor

template<Integer S, Integer N, Integer M>
Real Sandals::RungeKutta< S, N, M >::m_safety_factor {0.9}
private

Safety factor for adaptive step \( f \).

◆ m_system

template<Integer S, Integer N, Integer M>
System Sandals::RungeKutta< S, N, M >::m_system
private

ODE/DAE system object pointer.

◆ m_tableau

template<Integer S, Integer N, Integer M>
Tableau<S> Sandals::RungeKutta< S, N, M >::m_tableau
private

Butcher tableau of the Runge-Kutta method.

◆ m_verbose

template<Integer S, Integer N, Integer M>
bool Sandals::RungeKutta< S, N, M >::m_verbose {false}
private

Verbose mode boolean.

The documentation for this class was generated from the following file: