Sandals::RungeKutta< Real, S, N, M > Class Template Reference

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

#include <RungeKutta.hh>

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

Public Member Functions
	RungeKutta (const RungeKutta &)=delete
RungeKutta &	operator= (RungeKutta const &)=delete
	RungeKutta (Tableau< Real, S > const &t_tableau)
	RungeKutta (Tableau< Real, S > const &t_tableau, System t_system)
Type	type () const
bool	is_erk () const
bool	is_irk () const
bool	is_dirk () const
Tableau< Real, S > &	tableau ()
Tableau< Real, 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)
void	implicit_system (typename ImplicitWrapper< Real, N, M >::FunctionF F, typename ImplicitWrapper< Real, N, M >::FunctionJF JF_x, typename ImplicitWrapper< Real, N, M >::FunctionJF JF_x_dot, typename ImplicitWrapper< Real, N, M >::FunctionH h=ImplicitWrapper< Real, N, M >::DefaultH, typename ImplicitWrapper< Real, N, M >::FunctionJH Jh_x=ImplicitWrapper< Real, N, M >::DefaultJH, typename ImplicitWrapper< Real, N, M >::FunctionID in_domain=ImplicitWrapper< Real, N, M >::DefaultID)
void	implicit_system (std::string name, typename ImplicitWrapper< Real, N, M >::FunctionF F, typename ImplicitWrapper< Real, N, M >::FunctionJF JF_x, typename ImplicitWrapper< Real, N, M >::FunctionJF JF_x_dot, typename ImplicitWrapper< Real, N, M >::FunctionH h=ImplicitWrapper< Real, N, M >::DefaultH, typename ImplicitWrapper< Real, N, M >::FunctionJH Jh_x=ImplicitWrapper< Real, N, M >::DefaultJH, typename ImplicitWrapper< Real, N, M >::FunctionID in_domain=ImplicitWrapper< Real, N, M >::DefaultID)
void	explicit_system (typename ExplicitWrapper< Real, N, M >::FunctionF f, typename ExplicitWrapper< Real, N, M >::FunctionJF Jf_x, typename ExplicitWrapper< Real, N, M >::FunctionH h=ExplicitWrapper< Real, N, M >::DefaultH, typename ExplicitWrapper< Real, N, M >::FunctionJH Jh_x=ExplicitWrapper< Real, N, M >::DefaultJH, typename ExplicitWrapper< Real, N, M >::FunctionID in_domain=ExplicitWrapper< Real, N, M >::DefaultID)
void	explicit_system (std::string name, typename ExplicitWrapper< Real, N, M >::FunctionF f, typename ExplicitWrapper< Real, N, M >::FunctionJF Jf_x, typename ExplicitWrapper< Real, N, M >::FunctionH h=ExplicitWrapper< Real, N, M >::DefaultH, typename ExplicitWrapper< Real, N, M >::FunctionJH Jh_x=ExplicitWrapper< Real, N, M >::DefaultJH, typename ExplicitWrapper< Real, N, M >::FunctionID in_domain=ExplicitWrapper< Real, N, M >::DefaultID)
void	linear_system (typename LinearWrapper< Real, N, M >::FunctionE E, typename LinearWrapper< Real, N, M >::FunctionA A, typename LinearWrapper< Real, N, M >::FunctionB b, typename LinearWrapper< Real, N, M >::FunctionH h=LinearWrapper< Real, N, M >::DefaultH, typename LinearWrapper< Real, N, M >::FunctionJH Jh_x=LinearWrapper< Real, N, M >::DefaultJH, typename LinearWrapper< Real, N, M >::FunctionID in_domain=LinearWrapper< Real, N, M >::DefaultID)
void	linear_system (std::string name, typename LinearWrapper< Real, N, M >::FunctionE E, typename LinearWrapper< Real, N, M >::FunctionA A, typename LinearWrapper< Real, N, M >::FunctionB b, typename LinearWrapper< Real, N, M >::FunctionH h=LinearWrapper< Real, N, M >::DefaultH, typename LinearWrapper< Real, N, M >::FunctionJH Jh_x=LinearWrapper< Real, N, M >::DefaultJH, typename LinearWrapper< Real, N, M >::FunctionID in_domain=LinearWrapper< Real, N, M >::DefaultID)
void	semi_explicit_system (typename SemiExplicitWrapper< Real, N, M >::FunctionA A, typename SemiExplicitWrapper< Real, N, M >::FunctionTA TA_x, typename SemiExplicitWrapper< Real, N, M >::FunctionB b, typename SemiExplicitWrapper< Real, N, M >::FunctionJB Jb_x, typename SemiExplicitWrapper< Real, N, M >::FunctionH h=SemiExplicitWrapper< Real, N, M >::DefaultH, typename SemiExplicitWrapper< Real, N, M >::FunctionJH Jh_x=SemiExplicitWrapper< Real, N, M >::DefaultJH, typename SemiExplicitWrapper< Real, N, M >::FunctionID in_domain=SemiExplicitWrapper< Real, N, M >::DefaultID)
void	semi_explicit_system (std::string name, typename SemiExplicitWrapper< Real, N, M >::FunctionA A, typename SemiExplicitWrapper< Real, N, M >::FunctionTA TA_x, typename SemiExplicitWrapper< Real, N, M >::FunctionB b, typename SemiExplicitWrapper< Real, N, M >::FunctionJB Jb_x, typename SemiExplicitWrapper< Real, N, M >::FunctionH h=SemiExplicitWrapper< Real, N, M >::DefaultH, typename SemiExplicitWrapper< Real, N, M >::FunctionJH Jh_x=SemiExplicitWrapper< Real, N, M >::DefaultJH, typename SemiExplicitWrapper< Real, N, M >::FunctionID in_domain=SemiExplicitWrapper< Real, N, M >::DefaultID)
bool	has_system ()
Real	absolute_tolerance ()
void	absolute_tolerance (Real const t_absolute_tolerance)
Real	relative_tolerance ()
void	relative_tolerance (Real const t_relative_tolerance)
Real &	safety_factor ()
void	safety_factor (Real const t_safety_factor)
Real &	min_safety_factor ()
void	min_safety_factor (Real const t_min_safety_factor)
Real &	max_safety_factor ()
void	max_safety_factor (Real const t_max_safety_factor)
Real &	min_step ()
void	min_step (Real const t_min_step)
Integer &	max_substeps ()
void	max_substeps (Integer const 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 ()
FunctionSC	step_callback ()
void	step_callback (FunctionSC const &t_step_callback)
Real	projection_tolerance ()
void	projection_tolerance (Real const t_projection_tolerance)
Integer &	max_projection_iterations ()
void	max_projection_iterations (Integer const 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 const h_k) const
std::string	info () const
void	info (std::ostream &os)
bool	erk_explicit_step (VectorN const &x_old, Real const t_old, Real const h_old, VectorN &x_new, Real &h_new, MatrixK &K) const
void	erk_implicit_function (Integer const s, VectorN const &x, Real const t, Real const h, MatrixK const &K, VectorN &fun) const
void	erk_implicit_jacobian (Integer const s, VectorN const &x, Real const t, Real const h, MatrixK const &K, MatrixN &jac) const
bool	erk_implicit_step (VectorN const &x_old, Real const t_old, Real const h_old, VectorN &x_new, Real &h_new, MatrixK &K) const
void	irk_function (VectorN const &x, Real const t, Real const h, VectorK const &K, VectorK &fun) const
void	irk_jacobian (VectorN const &x, Real const t, Real const h, VectorK const &K, MatrixJ &jac) const
bool	irk_step (VectorN const &x_old, Real const t_old, Real const h_old, VectorN &x_new, Real &h_new, MatrixK &K) const
void	dirk_function (Integer n, VectorN const &x, Real const t, Real const h, MatrixK const &K, VectorN &fun) const
void	dirk_jacobian (Integer n, VectorN const &x, Real const t, Real const h, MatrixK const &K, MatrixN &jac) const
bool	dirk_step (VectorN const &x_old, Real const t_old, Real const h_old, VectorN &x_new, Real &h_new, MatrixK &K) const
bool	step (VectorN const &x_old, Real const t_old, Real const h_old, VectorN &x_new, Real &h_new, MatrixK &K) const
bool	advance (VectorN const &x_old, Real const t_old, Real const h_old, VectorN &x_new, Real &h_new) const
bool	solve (VectorX const &t_mesh, VectorN const &ics, Solution< Real, N, M > &sol) const
bool	adaptive_solve (VectorX const &t_mesh, VectorN const &ics, Solution< Real, N, M > &sol) const
bool	project (VectorN const &x, Real const t, VectorN &x_projected) const
bool	project_ics (VectorN const &x, Real const 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) const

Public Attributes
const Real	SQRT_EPSILON {std::sqrt(EPSILON)}

Private Types
using	VectorX = Eigen::Vector<Real, Eigen::Dynamic>
using	MatrixX = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>
using	VectorK = Eigen::Vector<Real, N*S>
using	MatrixK = Eigen::Matrix<Real, N, S>
using	MatrixJK = Eigen::Matrix<Real, N*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<Real, N>
using	NewtonK = Optimist::RootFinder::Newton<Real, N*S>
using	VectorS = typename Tableau<Real, S>::Vector
using	MatrixS = typename Tableau<Real, S>::Matrix
using	VectorN = typename Implicit<Real, N, M>::VectorF
using	MatrixN = typename Implicit<Real, N, M>::MatrixJF
using	VectorM = typename Implicit<Real, N, M>::VectorH
using	MatrixM = typename Implicit<Real, N, M>::MatrixJH
using	FunctionSC = std::function<void(Integer const, VectorX const &, Real const)>

Private Attributes
NewtonX	m_newtonX
NewtonK	m_newtonK
Eigen::FullPivLU< MatrixP >	m_lu
Tableau< Real, S >	m_tableau
System	m_system
Real	m_absolute_tolerance {1e-6}
Real	m_relative_tolerance {1e-3}
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}
FunctionSC	m_step_callback {nullptr}
Real	m_projection_tolerance {EPSILON_HIGH}
Integer	m_max_projection_iterations {5}
bool	m_projection {true}

Detailed Description

template<typename Real, Integer S, Integer N, Integer M>
class Sandals::RungeKutta< Real, 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) \\ \mathbf{K}_3 = h_k \mathbf{f} \left(\mathbf{x}_k + a_{31} \mathbf{K}_1 + a_{32} \mathbf{K}_2, t_k + h_k c_3\right) \\[-0.5em] \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{.} \]

In matrix form, this can be expressed as

\[ \begin{bmatrix} \mathbf{I} & \mathbf{0} & \cdots & \cdots & \mathbf{0} \\ a_{21} & \mathbf{I} & \ddots & \mathbf{0} & \vdots \\ a_{31} & a_{32} & \mathbf{I} & \ddots & \vdots \\ \vdots & \vdots & \ddots & \ddots & \mathbf{0} \\ a_{s1} & a_{s2} & \cdots & a_{s,s-1} & \mathbf{I} \end{bmatrix} \begin{bmatrix} \mathbf{K}_1 \\ \mathbf{K}_2 \\ \mathbf{K}_3 \\ \vdots \\ \mathbf{K}_s \end{bmatrix} = h_k \begin{bmatrix} \mathbf{f} \left(\mathbf{x}_k, t_k\right) \\ \mathbf{f} \left(\mathbf{x}_k + a_{21} \mathbf{K}_1, t_k + h_k c_2\right) \\ \mathbf{f} \left(\mathbf{x}_k + a_{31} \mathbf{K}_1 + a_{32} \mathbf{K}_2, t_k + h_k c_3\right) \\ \vdots \\ \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{bmatrix} \text{.} \]

Derivatives of \(\mathbf{K}\) with respect to \(\mathbf{x}_k\)

For explicit systems, the derivatives of the intermediate variables \(\mathbf{K}\) with respect to the states \(\mathbf{x}\) can be computed by the differentiation of the stage equations, which yields

\[ \begin{array}{l} \begin{cases} \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k} = h_k \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}(\mathbf{x}_k, t_k) \\ \displaystyle\frac{\partial\mathbf{K}_2}{\partial\mathbf{x}_k} = h_k \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + a_{21} \mathbf{K}_1, t_k + h_k c_2\right) \left(\mathbf{I} + a_{21} \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k}\right) \\ \displaystyle\frac{\partial\mathbf{K}_3}{\partial\mathbf{x}_k} = h_k \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + a_{31} \mathbf{K}_1 + a_{32} \mathbf{K}_2, t_k + h_k c_3\right) \left(\mathbf{I} + a_{31} \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k} + a_{32} \displaystyle\frac{\partial\mathbf{K}_2}{\partial\mathbf{x}_k}\right) \\[-0.5em] \vdots \\[-0.5em] \displaystyle\frac{\partial\mathbf{K}_s}{\partial\mathbf{x}_k} = h_k \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s-1} a_{sj} \mathbf{K}_j, t_k + h_k c_s\right) \left(\mathbf{I} + \displaystyle\sum_{j=1}^{s-1} a_{sj} \displaystyle\frac{\partial\mathbf{K}_j}{\partial\mathbf{x}_k}\right) \end{cases} \end{array} \text{.} \]

In matrix form, this can be expressed as

\[ \begin{bmatrix} \mathbf{I} & \mathbf{0} & \cdots & \cdots & \mathbf{0} \\ -h_k a_{21} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & \mathbf{I} & \ddots & \mathbf{0} & \vdots \\ -h_k a_{31} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & -h_k a_{32} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & \mathbf{I} & \ddots & \vdots \\ \vdots & \vdots & \ddots & \ddots & \mathbf{0} \\ -h_k a_{s1} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & -h_k a_{s2} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & \cdots & -h_k a_{s,s-1} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & \mathbf{I} \end{bmatrix} \begin{bmatrix} \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k} \\ \displaystyle\frac{\partial\mathbf{K}_2}{\partial\mathbf{x}_k} \\ \displaystyle\frac{\partial\mathbf{K}_3}{\partial\mathbf{x}_k} \\ \vdots \\ \displaystyle\frac{\partial\mathbf{K}_s}{\partial\mathbf{x}_k} \end{bmatrix} = h_k \begin{bmatrix} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}(\mathbf{x}_k, t_k) \\ \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + a_{21} \mathbf{K}_1, t_k + h_k c_2\right) \\ \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + a_{31} \mathbf{K}_1 + a_{32} \mathbf{K}_2, t_k + h_k c_3\right) \\ \vdots\\ \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s-1} a_{sj} \mathbf{K}_j, t_k + h_k c_s\right) \end{bmatrix} \text{.} \]

Derivatives of \(\mathbf{x}_{k+1}\) with respect to \(\mathbf{x}_k\)

The derivatives of the next state \(\mathbf{x}_{k+1}\) with respect to the states \(\mathbf{x}\) can be computed as

\[ \displaystyle\frac{\partial\mathbf{x}_{k+1}}{\partial\mathbf{x}_k} = \mathbf{I} + h_k \displaystyle\sum_{i=1}^s b_i \frac{\partial\mathbf{K}_i}{\partial\mathbf{x}_k} \text{.} \]

where the derivatives of the intermediate variables \(\mathbf{K}_i\) with respect to the states \(\mathbf{x}_k\) are computed as described above.

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, \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, \displaystyle\frac{\mathbf{K}_2}{h_k}, t_k + h_k c_2\right) = \mathbf{0} \\ \mathbf{F}_3 \left(\mathbf{x}_k + a_{31} \mathbf{K}_1 + a_{32} \mathbf{K}_2, \displaystyle\frac{\mathbf{K}_3}{h_k}, t_k + h_k c_3\right) = \mathbf{0} \\[-0.5em] \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.

Derivatives of \(\mathbf{K}\) with respect to \(\mathbf{x}_k\)

For implicit systems, the derivatives of the intermediate variables \(\mathbf{K}\) with respect to the states \(\mathbf{x}_k\) can be computed by differentiating the stage equations, which yields

\[ \begin{array}{l} \begin{cases} \displaystyle\frac{\partial\mathbf{F}_1}{\partial\mathbf{x}_k}\left(\mathbf{x}_k, \displaystyle\frac{\mathbf{K}_1}{h_k}, t_k + h_k c_1\right) + \frac{1}{h_k} \displaystyle\frac{\partial\mathbf{F}_1}{\partial\mathbf{K}_1}\left(\mathbf{x}_k, \displaystyle\frac{\mathbf{K}_1}{h_k}, t_k + h_k c_1\right) \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k} = \mathbf{0} \\ \displaystyle\frac{\partial\mathbf{F}_2}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + a_{21} \mathbf{K}_1, \displaystyle\frac{\mathbf{K}_2}{h_k}, t_k + h_k c_2\right) \left(\mathbf{I} + a_{21} \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k}\right) + \frac{1}{h_k} \displaystyle\frac{\partial\mathbf{F}_2}{\partial\mathbf{K}_2}\left(\mathbf{x}_k + a_{21} \mathbf{K}_1, \displaystyle\frac{\mathbf{K}_2}{h_k}, t_k + h_k c_2\right) \displaystyle\frac{\partial\mathbf{K}_2}{\partial\mathbf{x}_k} = \mathbf{0} \\ \displaystyle\frac{\partial\mathbf{F}_3}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + a_{31} \mathbf{K}_1 + a_{32} \mathbf{K}_2, \displaystyle\frac{\mathbf{K}_3}{h_k}, t_k + h_k c_3\right) \left(\mathbf{I} + a_{31} \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k} + a_{32} \displaystyle\frac{\partial\mathbf{K}_2}{\partial\mathbf{x}_k}\right) + \frac{1}{h_k} \displaystyle\frac{\partial\mathbf{F}_3}{\partial\mathbf{K}_3}\left(\mathbf{x}_k + a_{31} \mathbf{K}_1 + a_{32} \mathbf{K}_2, \displaystyle\frac{\mathbf{K}_3}{h_k}, t_k + h_k c_3\right) \displaystyle\frac{\partial\mathbf{K}_3}{\partial\mathbf{x}_k} = \mathbf{0} \\[-0.5em] \vdots \\[-0.5em] \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{x}_k}\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) \left(\mathbf{I} + \displaystyle\sum_{j=1}^{s-1} a_{sj} \displaystyle\frac{\partial\mathbf{K}_j}{\partial\mathbf{x}_k}\right) + \frac{1}{h_k} \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{K}_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) \displaystyle\frac{\partial\mathbf{K}_s}{\partial\mathbf{x}_k} = \mathbf{0} \end{cases} \end{array} \text{.} \]

In matrix form, this can be expressed as

\[ \begin{bmatrix} \displaystyle\frac{\partial\mathbf{F}_1}{\partial\mathbf{K}_1}\left(\cdot\right) & \mathbf{0} & \cdots & \cdots & \mathbf{0} \\ h_k a_{21} \displaystyle\frac{\partial\mathbf{F}_2}{\partial\mathbf{x}_k}\left(\cdot\right) & \displaystyle\frac{\partial\mathbf{F}_2}{\partial\mathbf{K}_2}\left(\cdot\right) & \ddots & \mathbf{0} & \vdots \\ h_k a_{31} \displaystyle\frac{\partial\mathbf{F}_3}{\partial\mathbf{x}_k}\left(\cdot\right) & h_k a_{32} \displaystyle\frac{\partial\mathbf{F}_3}{\partial\mathbf{x}_k}\left(\cdot\right) & \displaystyle\frac{\partial\mathbf{F}_3}{\partial\mathbf{K}_3}\left(\cdot\right) & \ddots & \vdots \\ \vdots & \vdots & \ddots & \ddots & \mathbf{0} \\ h_k a_{s1} \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{x}_k}\left(\cdot\right) & h_k a_{s2} \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{x}_k}\left(\cdot\right) & \cdots & h_k a_{s,s-1} \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{x}_k}\left(\cdot\right) & \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{K}_s}\left(\cdot\right) \end{bmatrix} \begin{bmatrix} \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k} \\ \displaystyle\frac{\partial\mathbf{K}_2}{\partial\mathbf{x}_k} \\ \displaystyle\frac{\partial\mathbf{K}_3}{\partial\mathbf{x}_k} \\ \vdots \\ \displaystyle\frac{\partial\mathbf{K}_s}{\partial\mathbf{x}_k} \end{bmatrix} = -h_k \begin{bmatrix} \displaystyle\frac{\partial\mathbf{F}_1}{\partial\mathbf{x}_k}\left(\mathbf{x}_k, \displaystyle\frac{\mathbf{K}_1}{h_k}, t_k + h_k c_1\right) \\ \displaystyle\frac{\partial\mathbf{F}_2}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + a_{21} \mathbf{K}_1, \displaystyle\frac{\mathbf{K}_2}{h_k}, t_k + h_k c_2\right) \\ \displaystyle\frac{\partial\mathbf{F}_3}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + a_{31} \mathbf{K}_1 + a_{32} \mathbf{K}_2, \displaystyle\frac{\mathbf{K}_3}{h_k}, t_k + h_k c_3\right) \\ \vdots\\ \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{x}_k}\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) \end{bmatrix} \text{.} \]

Derivatives of \(\mathbf{x}_{k+1}\) with respect to \(\mathbf{x}_k\)

NOTE: The derivatives of the next state \(\mathbf{x}_{k+1}\) with respect to the states \(\mathbf{x}_k\) can be computed as previously described for explicit systems, as the structure of the ERK method remains the same.

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) \\ \mathbf{K}_3 = h_k \mathbf{f} \left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{3j} \mathbf{K}_j, t_k + h_k c_3\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{.} \]

Derivatives of \(\mathbf{K}\) with respect to \(\mathbf{x}_k\)

For explicit systems, the derivatives of the intermediate variables \(\mathbf{K}\) with respect to the states \(\mathbf{x}_k\) can be computed by differentiating the stage equations, which yields

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

In matrix form, this can be expressed as

\[ \begin{bmatrix} \mathbf{I} - h_k a_{11} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & -h_k a_{12} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & -h_k a_{13} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & \cdots & -h_k a_{1s} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) \\ -h_k a_{21} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & \mathbf{I} - h_k a_{22} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & -h_k a_{23} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & \cdots & \vdots \\ -h_k a_{31} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & -h_k a_{32} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & \mathbf{I} - h_k a_{33} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & \ddots & \vdots \\ \vdots & \vdots & \ddots & \ddots & -h_k a_{s,s-1} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) \\ -h_k a_{s1} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & -h_k a_{s2} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & \cdots & -h_k a_{s,s-1} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & \mathbf{I} - h_k a_{ss} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) \end{bmatrix} \begin{bmatrix} \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k} \\ \displaystyle\frac{\partial\mathbf{K}_2}{\partial\mathbf{x}_k} \\ \displaystyle\frac{\partial\mathbf{K}_3}{\partial\mathbf{x}_k} \\ \vdots \\ \displaystyle\frac{\partial\mathbf{K}_s}{\partial\mathbf{x}_k} \end{bmatrix} = h_k \begin{bmatrix} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{1j} \mathbf{K}_j, t_k + h_k c_1\right) \\ \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{2j} \mathbf{K}_j, t_k + h_k c_2\right) \\ \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{3j} \mathbf{K}_j, t_k + h_k c_3\right) \\ \vdots \\ \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{sj} \mathbf{K}_j, t_k + h_k c_s\right) \end{bmatrix} \text{.} \]

Derivatives of \(\mathbf{x}_{k+1}\) with respect to \(\mathbf{x}_k\)

The derivatives of the next state \(\mathbf{x}_{k+1}\) with respect to the states \(\mathbf{x}\) can be computed as

\[ \displaystyle\frac{\partial\mathbf{x}_{k+1}}{\partial\mathbf{x}_k} = \mathbf{I} + h_k \displaystyle\sum_{i=1}^s b_i \frac{\partial\mathbf{K}_i}{\partial\mathbf{x}_k} \text{.} \]

where the derivatives of the intermediate variables \(\mathbf{K}_i\) with respect to the states \(\mathbf{x}_k\) are computed as described above.

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} \\ \mathbf{F}_3 \left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{3j} \mathbf{K}_j, \displaystyle\frac{\mathbf{K}_3}{h_k}, t_k + h_k c_3\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.

Derivatives of \(\mathbf{K}\) with respect to \(\mathbf{x}_k\)

For implicit systems, the derivatives of the intermediate variables \(\mathbf{K}\) with respect to the states \(\mathbf{x}_k\) can be computed by differentiating the stage equations, which yields

\[ \begin{array}{l} \begin{cases} \displaystyle\frac{\partial\mathbf{F}_1}{\partial\mathbf{x}_k}\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) \left(\mathbf{I} + \displaystyle\sum_{j=1}^{s} a_{1j} \displaystyle\frac{\partial\mathbf{K}_j}{\partial\mathbf{x}_k}\right) + \frac{1}{h_k} \displaystyle\frac{\partial\mathbf{F}_1}{\partial\mathbf{K}_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) \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k} = \mathbf{0} \\ \displaystyle\frac{\partial\mathbf{F}_2}{\partial\mathbf{x}_k}\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) \left(\mathbf{I} + \displaystyle\sum_{j=1}^{s} a_{2j} \displaystyle\frac{\partial\mathbf{K}_j}{\partial\mathbf{x}_k}\right) + \frac{1}{h_k} \displaystyle\frac{\partial\mathbf{F}_2}{\partial\mathbf{K}_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) \displaystyle\frac{\partial\mathbf{K}_2}{\partial\mathbf{x}_k} = \mathbf{0} \\ \displaystyle\frac{\partial\mathbf{F}_3}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{3j} \mathbf{K}_j, \displaystyle\frac{\mathbf{K}_3}{h_k}, t_k + h_k c_3\right) \left(\mathbf{I} + \displaystyle\sum_{j=1}^{s} a_{3j} \displaystyle\frac{\partial\mathbf{K}_j}{\partial\mathbf{x}_k}\right) + \frac{1}{h_k} \displaystyle\frac{\partial\mathbf{F}_3}{\partial\mathbf{K}_3}\left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{3j} \mathbf{K}_j, \displaystyle\frac{\mathbf{K}_3}{h_k}, t_k + h_k c_3\right) \displaystyle\frac{\partial\mathbf{K}_3}{\partial\mathbf{x}_k} = \mathbf{0} \\[-0.5em] \vdots \\[-0.5em] \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{x}_k}\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) \left(\mathbf{I} + \displaystyle\sum_{j=1}^{s} a_{sj} \displaystyle\frac{\partial\mathbf{K}_j}{\partial\mathbf{x}_k}\right) + \frac{1}{h_k} \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{K}_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) \displaystyle\frac{\partial\mathbf{K}_s}{\partial\mathbf{x}_k} = \mathbf{0} \end{cases} \end{array} \text{.} \]

In matrix form, this can be expressed as

\[ \begin{bmatrix} \displaystyle\frac{\partial\mathbf{F}_1}{\partial\mathbf{K}_1}\left(\cdot\right) & h_k a_{12} \displaystyle\frac{\partial\mathbf{F}_1}{\partial\mathbf{x}_k}\left(\cdot\right) & h_k a_{13} \displaystyle\frac{\partial\mathbf{F}_1}{\partial\mathbf{x}_k}\left(\cdot\right) & \cdots & h_k a_{1s} \displaystyle\frac{\partial\mathbf{F}_1}{\partial\mathbf{x}_k}\left(\cdot\right) \\ h_k a_{21} \displaystyle\frac{\partial\mathbf{F}_2}{\partial\mathbf{x}_k}\left(\cdot\right) & \displaystyle\frac{\partial\mathbf{F}_2}{\partial\mathbf{K}_2}\left(\cdot\right) & h_k a_{23} \displaystyle\frac{\partial\mathbf{F}_2}{\partial\mathbf{x}_k}\left(\cdot\right) & \cdots & \vdots \\ h_k a_{31} \displaystyle\frac{\partial\mathbf{F}_3}{\partial\mathbf{x}_k}\left(\cdot\right) & h_k a_{32} \displaystyle\frac{\partial\mathbf{F}_3}{\partial\mathbf{x}_k}\left(\cdot\right) & \displaystyle\frac{\partial\mathbf{F}_3}{\partial\mathbf{K}_3}\left(\cdot\right) & \ddots & \vdots \\ \vdots & \vdots & \ddots & \ddots & h_k a_{s,s-1} \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{x}_k}\left(\cdot\right) \\ h_k a_{s1} \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{x}_k}\left(\cdot\right) & h_k a_{s2} \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{x}_k}\left(\cdot\right) & \cdots & h_k a_{s,s-1} \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{x}_k}\left(\cdot\right) & \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{K}_s}\left(\cdot\right) \end{bmatrix} \begin{bmatrix} \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k} \\ \displaystyle\frac{\partial\mathbf{K}_2}{\partial\mathbf{x}_k} \\ \displaystyle\frac{\partial\mathbf{K}_3}{\partial\mathbf{x}_k} \\ \vdots \\ \displaystyle\frac{\partial\mathbf{K}_s}{\partial\mathbf{x}_k} \end{bmatrix} = -h_k \begin{bmatrix} \displaystyle\frac{\partial\mathbf{F}_1}{\partial\mathbf{x}_k}\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) \\ \displaystyle\frac{\partial\mathbf{F}_2}{\partial\mathbf{x}_k}\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) \\ \displaystyle\frac{\partial\mathbf{F}_3}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{3j} \mathbf{K}_j, \displaystyle\frac{\mathbf{K}_3}{h_k}, t_k + h_k c_3\right) \\ \vdots \\ \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{x}_k}\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) \end{bmatrix} \text{.} \]

Derivatives of \(\mathbf{x}_{k+1}\) with respect to \(\mathbf{x}_k\)

NOTE: The derivatives of the next state \(\mathbf{x}_{k+1}\) with respect to the states \(\mathbf{x}_k\) can be computed as previously described for explicit systems, as the structure of the IRK method remains the same.

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{.} \]

Derivatives of \(\mathbf{K}\) with respect to \(\mathbf{x}_k\)

For explicit systems, the derivatives of the intermediate variables \(\mathbf{K}\) with respect to the states \(\mathbf{x}_k\) can be computed by differentiating the stage equations, which yields

\[ \begin{array}{l} \begin{cases} \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k} = h_k \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + a_{11} \mathbf{K}_1, t_k + h_k c_1\right) \left(\mathbf{I} + a_{11} \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k}\right) \\ \displaystyle\frac{\partial\mathbf{K}_2}{\partial\mathbf{x}_k} = h_k \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + a_{21} \mathbf{K}_1 + a_{22} \mathbf{K}_2, t_k + h_k c_2\right) \left(\mathbf{I} + a_{21} \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k} + a_{22} \displaystyle\frac{\partial\mathbf{K}_2}{\partial\mathbf{x}_k}\right) \\ \displaystyle\frac{\partial\mathbf{K}_3}{\partial\mathbf{x}_k} = h_k \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + a_{31} \mathbf{K}_1 + a_{32} \mathbf{K}_2 + a_{33} \mathbf{K}_3, t_k + h_k c_3\right) \left(\mathbf{I} + a_{31} \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k} + a_{32} \displaystyle\frac{\partial\mathbf{K}_2}{\partial\mathbf{x}_k} + a_{33} \displaystyle\frac{\partial\mathbf{K}_3}{\partial\mathbf{x}_k}\right) \\[-0.5em] \vdots \\[-0.5em] \displaystyle\frac{\partial\mathbf{K}_s}{\partial\mathbf{x}_k} = h_k \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + \displaystyle\sum_{j=1}^{s} a_{sj} \mathbf{K}_j, t_k + h_k c_s\right) \left(\mathbf{I} + \displaystyle\sum_{j=1}^{s} a_{sj} \displaystyle\frac{\partial\mathbf{K}_j}{\partial\mathbf{x}_k}\right) \end{cases} \end{array} \text{.} \]

In matrix form, this can be expressed as

\[ \begin{bmatrix} \mathbf{I} - h_k a_{11} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} \\ -h_k a_{21} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & \mathbf{I} - h_k a_{22} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & \ddots & \mathbf{0} & \vdots \\ -h_k a_{31} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & -h_k a_{32} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & \mathbf{I} - h_k a_{33} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & \ddots & \vdots \\ \vdots & \vdots & \ddots & \ddots & \mathbf{0} \\ -h_k a_{s1} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & -h_k a_{s2} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & \cdots & -h_k a_{s,s-1} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) & \mathbf{I} - h_k a_{ss} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) \end{bmatrix} \begin{bmatrix} \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k} \\ \displaystyle\frac{\partial\mathbf{K}_2}{\partial\mathbf{x}_k} \\ \displaystyle\frac{\partial\mathbf{K}_3}{\partial\mathbf{x}_k} \\ \vdots \\ \displaystyle\frac{\partial\mathbf{K}_s}{\partial\mathbf{x}_k} \end{bmatrix} = h_k \begin{bmatrix} \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) \\ \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) \\ \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) \\ \vdots \\ \displaystyle\frac{\partial\mathbf{f}}{\partial\mathbf{x}_k}\left(\cdot\right) \end{bmatrix} \text{.} \]

Derivatives of \(\mathbf{x}_{k+1}\) with respect to \(\mathbf{x}_k\)

The derivatives of the next state \(\mathbf{x}_{k+1}\) with respect to the states \(\mathbf{x}\) can be computed as

\[ \displaystyle\frac{\partial\mathbf{x}_{k+1}}{\partial\mathbf{x}_k} = \mathbf{I} + h_k \displaystyle\sum_{i=1}^s b_i \frac{\partial\mathbf{K}_i}{\partial\mathbf{x}_k} \text{.} \]

where the derivatives of the intermediate variables \(\mathbf{K}_i\) with respect to the states \(\mathbf{x}_k\) are computed as described above.

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} \\ \mathbf{F}_3 \left(\mathbf{x}_k + a_{31} \mathbf{K}_1 + a_{32} \mathbf{K}_2 + a_{33} \mathbf{K}_3, \displaystyle\frac{\mathbf{K}_3}{h_k}, t_k + h_k c_3\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: 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.

Derivatives of \(\mathbf{K}\) with respect to \(\mathbf{x}_k\)

For implicit systems, the derivatives of the intermediate variables \(\mathbf{K}\) with respect to the states \(\mathbf{x}_k\) can be computed by differentiating the stage equations, which yields

\[ \begin{array}{l} \begin{cases} \displaystyle\frac{\partial\mathbf{F}_1}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + a_{11} \mathbf{K}_1, \displaystyle\frac{\mathbf{K}_1}{h_k}, t_k + h_k c_1\right) \left(\mathbf{I} + a_{11} \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k}\right) + \frac{1}{h_k} \displaystyle\frac{\partial\mathbf{F}_1}{\partial\mathbf{K}_1}\left(\mathbf{x}_k + a_{11} \mathbf{K}_1, \displaystyle\frac{\mathbf{K}_1}{h_k}, t_k + h_k c_1\right) \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k} = \mathbf{0} \\ \displaystyle\frac{\partial\mathbf{F}_2}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + a_{21} \mathbf{K}_1 + a_{22} \mathbf{K}_2, \displaystyle\frac{\mathbf{K}_2}{h_k}, t_k + h_k c_2\right) \left(\mathbf{I} + a_{21} \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k} + a_{22} \displaystyle\frac{\partial\mathbf{K}_2}{\partial\mathbf{x}_k}\right) + \frac{1}{h_k} \displaystyle\frac{\partial\mathbf{F}_2}{\partial\mathbf{K}_2}\left(\mathbf{x}_k + a_{21} \mathbf{K}_1 + a_{22} \mathbf{K}_2, \displaystyle\frac{\mathbf{K}_2}{h_k}, t_k + h_k c_2\right) \displaystyle\frac{\partial\mathbf{K}_2}{\partial\mathbf{x}_k} = \mathbf{0} \\[-0.5em] \displaystyle\frac{\partial\mathbf{F}_3}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + a_{31} \mathbf{K}_1 + a_{32} \mathbf{K}_2 + a_{33} \mathbf{K}_3, \displaystyle\frac{\mathbf{K}_3}{h_k}, t_k + h_k c_3\right) \left(\mathbf{I} + a_{31} \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k} + a_{32} \displaystyle\frac{\partial\mathbf{K}_2}{\partial\mathbf{x}_k} + a_{33} \displaystyle\frac{\partial\mathbf{K}_3}{\partial\mathbf{x}_k}\right) + \frac{1}{h_k} \displaystyle\frac{\partial\mathbf{F}_3}{\partial\mathbf{K}_3}\left(\mathbf{x}_k + a_{31} \mathbf{K}_1 + a_{32} \mathbf{K}_2 + a_{33} \mathbf{K}_3, \displaystyle\frac{\mathbf{K}_3}{h_k}, t_k + h_k c_3\right) \displaystyle\frac{\partial\mathbf{K}_3}{\partial\mathbf{x}_k} = \mathbf{0} \\[-0.5em] \vdots \\[-0.5em] \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{x}_k}\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) \left(\mathbf{I} + \displaystyle\sum_{j=1}^{s} a_{sj} \displaystyle\frac{\partial\mathbf{K}_j}{\partial\mathbf{x}_k}\right) + \frac{1}{h_k} \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{K}_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) \displaystyle\frac{\partial\mathbf{K}_s}{\partial\mathbf{x}_k} = \mathbf{0} \end{cases} \end{array} \text{.} \]

In matrix form, this can be expressed as

\[ \begin{bmatrix} h_k a_{11} \displaystyle\frac{\partial\mathbf{F}_1}{\partial\mathbf{x}_k}\left(\cdot\right) & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} \\ -h_k a_{21} \displaystyle\frac{\partial\mathbf{F}_2}{\partial\mathbf{x}_k}\left(\cdot\right) & h_k a_{22} \displaystyle\frac{\partial\mathbf{F}_2}{\partial\mathbf{x}_k}\left(\cdot\right) & \ddots & \mathbf{0} & \vdots \\ -h_k a_{31} \displaystyle\frac{\partial\mathbf{F}_3}{\partial\mathbf{x}_k}\left(\cdot\right) & -h_k a_{32} \displaystyle\frac{\partial\mathbf{F}_3}{\partial\mathbf{x}_k}\left(\cdot\right) & h_k a_{33} \displaystyle\frac{\partial\mathbf{F}_3}{\partial\mathbf{x}_k}\left(\cdot\right) & \ddots & \vdots \\ \vdots & \vdots & \ddots & \ddots & \mathbf{0} \\ -h_k a_{s1} \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{x}_k}\left(\cdot\right) & -h_k a_{s2} \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{x}_k}\left(\cdot\right) & \cdots & -h_k a_{s,s-1} \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{x}_k}\left(\cdot\right) & \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{K}_s}\left(\cdot\right) \end{bmatrix} \begin{bmatrix} \displaystyle\frac{\partial\mathbf{K}_1}{\partial\mathbf{x}_k} \\ \displaystyle\frac{\partial\mathbf{K}_2}{\partial\mathbf{x}_k} \\ \displaystyle\frac{\partial\mathbf{K}_3}{\partial\mathbf{x}_k} \\ \vdots \\ \displaystyle\frac{\partial\mathbf{K}_s}{\partial\mathbf{x}_k} \end{bmatrix} = -h_k \begin{bmatrix} \displaystyle\frac{\partial\mathbf{F}_1}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + a_{11} \mathbf{K}_1, \displaystyle\frac{\mathbf{K}_1}{h_k}, t_k + h_k c_1\right) \\ \displaystyle\frac{\partial\mathbf{F}_2}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + a_{21} \mathbf{K}_1 + a_{22} \mathbf{K}_2, \displaystyle\frac{\mathbf{K}_2}{h_k}, t_k + h_k c_2\right) \\ \displaystyle\frac{\partial\mathbf{F}_3}{\partial\mathbf{x}_k}\left(\mathbf{x}_k + a_{31} \mathbf{K}_1 + a_{32} \mathbf{K}_2 + a_{33} \mathbf{K}_3, \displaystyle\frac{\mathbf{K}_3}{h_k}, t_k + h_k c_3\right) \\ \vdots \\ \displaystyle\frac{\partial\mathbf{F}_s}{\partial\mathbf{x}_k}\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) \end{bmatrix} \text{.} \]

Derivatives of \(\mathbf{x}_{k+1}\) with respect to \(\mathbf{x}_k\)

NOTE: The derivatives of the next state \(\mathbf{x}_{k+1}\) with respect to the states \(\mathbf{x}_k\) can be computed as previously described for explicit systems, as the structure of the IRK method remains the same.

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

Real	The scalar number type.
S	The number of stages of the Runge-Kutta method.
N	The dimension of the ODE/DAE system.
M	The dimension of the invariants manifold.

Member Typedef Documentation

FunctionSC

template<typename Real, Integer S, Integer N, Integer M>

using Sandals::RungeKutta< Real, S, N, M >::FunctionSC = std::function<void(Integer const, VectorX const &, Real const)>

private

Step callback function type.

MatrixJ

template<typename Real, Integer S, Integer N, Integer M>

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

private

Templetized matrix type.

MatrixJK

template<typename Real, Integer S, Integer N, Integer M>

using Sandals::RungeKutta< Real, S, N, M >::MatrixJK = Eigen::Matrix<Real, N*N, S>

private

Templetized matrix type.

MatrixK

template<typename Real, Integer S, Integer N, Integer M>

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

private

Templetized matrix type.

MatrixM

template<typename Real, Integer S, Integer N, Integer M>

using Sandals::RungeKutta< Real, S, N, M >::MatrixM = typename Implicit<Real, N, M>::MatrixJH

private

Templetized matrix type.

MatrixN

template<typename Real, Integer S, Integer N, Integer M>

using Sandals::RungeKutta< Real, S, N, M >::MatrixN = typename Implicit<Real, N, M>::MatrixJF

private

Templetized matrix type.

MatrixP

template<typename Real, Integer S, Integer N, Integer M>

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

private

Templetized matrix type.

MatrixS

template<typename Real, Integer S, Integer N, Integer M>

using Sandals::RungeKutta< Real, S, N, M >::MatrixS = typename Tableau<Real, S>::Matrix

private

Templetized matrix type.

MatrixX

template<typename Real, Integer S, Integer N, Integer M>

using Sandals::RungeKutta< Real, S, N, M >::MatrixX = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>

private

\( N \times N \) matrix of Real number type.

NewtonK

template<typename Real, Integer S, Integer N, Integer M>

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

private

Templetized Newton solver for IRK methods.

NewtonX

template<typename Real, Integer S, Integer N, Integer M>

using Sandals::RungeKutta< Real, S, N, M >::NewtonX = Optimist::RootFinder::Newton<Real, N>

private

Templetized Newton solver for ERK and DIRK methods.

System

template<typename Real, Integer S, Integer N, Integer M>

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

Shared pointer to an implicit ODE/DAE system.

Time

template<typename Real, Integer S, Integer N, Integer M>

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

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

Type

template<typename Real, Integer S, Integer N, Integer M>

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

Runge-Kutta type enumeration.

VectorK

template<typename Real, Integer S, Integer N, Integer M>

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

private

Templetized vector type.

VectorM

template<typename Real, Integer S, Integer N, Integer M>

using Sandals::RungeKutta< Real, S, N, M >::VectorM = typename Implicit<Real, N, M>::VectorH

private

Templetized vector type.

VectorN

template<typename Real, Integer S, Integer N, Integer M>

using Sandals::RungeKutta< Real, S, N, M >::VectorN = typename Implicit<Real, N, M>::VectorF

private

Templetized vector type.

VectorP

template<typename Real, Integer S, Integer N, Integer M>

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

private

Templetized vector type.

VectorS

template<typename Real, Integer S, Integer N, Integer M>

using Sandals::RungeKutta< Real, S, N, M >::VectorS = typename Tableau<Real, S>::Vector

private

Templetized vector type.

VectorX

template<typename Real, Integer S, Integer N, Integer M>

using Sandals::RungeKutta< Real, S, N, M >::VectorX = Eigen::Vector<Real, Eigen::Dynamic>

private

\( N \times 1 \) vector of Real number type (column vector).

Constructor & Destructor Documentation

RungeKutta() [1/3]

template<typename Real, Integer S, Integer N, Integer M>

Sandals::RungeKutta< Real, S, N, M >::RungeKutta ( const RungeKutta< Real, S, N, M > & )

delete

Copy constructor for the timer.

RungeKutta() [2/3]

template<typename Real, Integer S, Integer N, Integer M>

Sandals::RungeKutta< Real, S, N, M >::RungeKutta ( Tableau< Real, S > const & t_tableau )

inline

Class constructor for the Runge-Kutta method.

Parameters

[in]

t_tableau

The Tableau reference.

RungeKutta() [3/3]

template<typename Real, Integer S, Integer N, Integer M>

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

inline

Class constructor for the Runge-Kutta method.

Parameters

[in]	t_tableau	The Tableau reference.
[in]	t_system	The ODE/DAE system shared pointer.

Member Function Documentation

A()

template<typename Real, Integer S, Integer N, Integer M>

MatrixS Sandals::RungeKutta< Real, 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<typename Real, Integer S, Integer N, Integer M>

Real Sandals::RungeKutta< Real, 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<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::absolute_tolerance ( Real const t_absolute_tolerance )

inline

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

Parameters

[in]

t_absolute_tolerance

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

adaptive()

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::adaptive ( bool t_adaptive )

inline

Set the adaptive step mode.

Parameters

[in]

t_adaptive

The adaptive step mode.

adaptive_mode()

template<typename Real, Integer S, Integer N, Integer M>

bool Sandals::RungeKutta< Real, S, N, M >::adaptive_mode ( )

inline

Get the adaptive step mode.

Returns

The adaptive step mode.

adaptive_solve()

template<typename Real, Integer S, Integer N, Integer M>

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

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_mesh	Independent variable (or time) mesh \( \mathbf{t} \).
[in]	ics	Initial conditions \( \mathbf{x}(t = 0) \).
[out]	sol	The solution of the system over the mesh of independent variable.

Returns

True if the system is successfully solved, false otherwise.

Warning

Do not use the solution for internal backtracking, as the step callback may directly modify the solution.

advance()

template<typename Real, Integer S, Integer N, Integer M>

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

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_old	States \( \mathbf{x}_k \) at the \( k \)-th step.
[in]	t_old	Independent variable (or time) \( t_k \) at the \( k \)-th step.
[in]	h_old	Advancing step \( h_k \) at the \( k \)-th step.
[out]	x_new	Computed states \( \mathbf{x}_{k+1} \) at the \( (k+1) \)-th step.
[out]	h_new	The suggested step \( h_{k+1}^\star \) for the next advancing step.

Returns

True if the step is successfully computed, false otherwise.

b()

template<typename Real, Integer S, Integer N, Integer M>

VectorS Sandals::RungeKutta< Real, 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<typename Real, Integer S, Integer N, Integer M>

VectorS Sandals::RungeKutta< Real, 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<typename Real, Integer S, Integer N, Integer M>

VectorS Sandals::RungeKutta< Real, 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<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::dirk_function ( Integer n, VectorN const & x, Real const t, Real const 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]	n	Equation index \( n \).
[in]	x	States \( \mathbf{x} \).
[in]	t	Independent variable (or time) \( t \).
[in]	h	Advancing step \( h \).
[in]	K	Variables \( \tilde{\mathbf{K}} = h \mathbf{K} \) of the system to be solved.
[out]	fun	The residual of system to be solved.

dirk_jacobian()

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::dirk_jacobian ( Integer n, VectorN const & x, Real const t, Real const 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]	n	Equation index \( n \).
[in]	x	States \( \mathbf{x} \).
[in]	t	Independent variable (or time) \( t \).
[in]	h	Advancing step \( h \).
[in]	K	Variables \( h \mathbf{K} \) of the system to be solved.
[out]	jac	The Jacobian of system to be solved.

dirk_step()

template<typename Real, Integer S, Integer N, Integer M>

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

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_old	States \( \mathbf{x}_k \) at the \( k \)-th step.
[in]	t_old	Independent variable (or time) \( t_k \) at the \( k \)-th step.
[in]	h_old	Advancing step \( h_k \) at the \( k \)-th step.
[out]	x_new	Computed states \( \mathbf{x}_{k+1} \) at the \( (k+1) \)-th step.
[out]	h_new	The suggested step \( h_{k+1}^\star \) for the next advancing step.
[out]	K	The \( \mathbf{K} \) variables of the Runge-Kutta method.

Returns

True if the step is successfully computed, false otherwise.

disable_adaptive_mode()

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::disable_adaptive_mode ( )

inline

Disable the adaptive step mode.

disable_projection()

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::disable_projection ( )

inline

Disable the projection mode.

disable_reverse_mode()

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::disable_reverse_mode ( )

inline

Disable the time reverse mode.

disable_verbose_mode()

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::disable_verbose_mode ( )

inline

Disable the verbose mode.

enable_adaptive_mode()

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::enable_adaptive_mode ( )

inline

Enable the adaptive step mode.

enable_projection()

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::enable_projection ( )

inline

Enable the projection mode.

enable_reverse_mode()

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::enable_reverse_mode ( )

inline

Enable the time reverse mode.

enable_verbose_mode()

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::enable_verbose_mode ( )

inline

Enable the verbose mode.

erk_explicit_step()

template<typename Real, Integer S, Integer N, Integer M>

bool Sandals::RungeKutta< Real, S, N, M >::erk_explicit_step ( VectorN const & x_old, Real const t_old, Real const h_old, VectorN & x_new, Real & h_new, MatrixK & K ) 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_old	States \( \mathbf{x}_k \) at the \( k \)-th step.
[in]	t_old	Independent variable (or time) \( t_k \) at the \( k \)-th step.
[in]	h_old	Advancing step \( h_k \) at the \( k \)-th step.
[out]	x_new	Computed states \( \mathbf{x}_{k+1} \) at the \( (k+1) \)-th step.
[out]	h_new	The suggested step \( h_{k+1}^\star \) for the next advancing step.
[out]	K	The \( \mathbf{K} \) variables of the Runge-Kutta method.

Returns

True if the step is successfully computed, false otherwise.

erk_implicit_function()

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::erk_implicit_function ( Integer const s, VectorN const & x, Real const t, Real const 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]	s	Stage index \( s \).
[in]	x	States \( \mathbf{x} \).
[in]	t	Independent variable (or time) \( t \).
[in]	h	Advancing step \( h \).
[in]	K	Variables \( \tilde{\mathbf{K}} = h \mathbf{K} \) of the system to be solved.
[out]	fun	The residual of system to be solved.

erk_implicit_jacobian()

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::erk_implicit_jacobian ( Integer const s, VectorN const & x, Real const t, Real const 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]	s	Stage index \( s \).
[in]	x	States \( \mathbf{x} \).
[in]	t	Independent variable (or time) \( t \).
[in]	h	Advancing step \( h \).
[in]	K	Variables \( h \mathbf{K} \) of the system to be solved.
[out]	jac	The Jacobian of system to be solved.

erk_implicit_step()

template<typename Real, Integer S, Integer N, Integer M>

bool Sandals::RungeKutta< Real, S, N, M >::erk_implicit_step ( VectorN const & x_old, Real const t_old, Real const h_old, VectorN & x_new, Real & h_new, MatrixK & K ) 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_old	States \( \mathbf{x}_k \) at the \( k \)-th step.
[in]	t_old	Independent variable (or time) \( t_k \) at the \( k \)-th step.
[in]	h_old	Advancing step \( h_k \) at the \( k \)-th step.
[out]	x_new	Computed states \( \mathbf{x}_{k+1} \) at the \( (k+1) \)-th step.
[out]	h_new	The suggested step \( h_{k+1}^\star \) for the next advancing step.
[out]	K	The \( \mathbf{K} \) variables of the Runge-Kutta method.

Returns

True if the step is successfully computed, false otherwise.

estimate_order()

template<typename Real, Integer S, Integer N, Integer M>

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

inline

Estimate the order of the Runge-Kutta method.

Parameters

[in]	t_mesh	The vector of time meshes with same initial and final time with fixed step.
[in]	ics	Initial conditions \( \mathbf{x}(t = 0) \).
[in]	sol	The analytical solution function.

Returns

The estimated order of the method.

estimate_step()

template<typename Real, Integer S, Integer N, Integer M>

Real Sandals::RungeKutta< Real, S, N, M >::estimate_step ( VectorN const & x, VectorN const & x_e, Real const 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]	x	States approximation \( \mathbf{x}_{k+1} \).
[in]	x_e	Embedded method's states approximation \( \hat{\mathbf{x}}_{k+1} \).
[in]	h_k	Actual advancing step \( h_k \).

Returns

The suggested step for the next integration step \( h_{k+1}^\star \).

explicit_system() [1/2]

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::explicit_system ( std::string name, typename ExplicitWrapper< Real, N, M >::FunctionF f, typename ExplicitWrapper< Real, N, M >::FunctionJF Jf_x, typename ExplicitWrapper< Real, N, M >::FunctionH h = ExplicitWrapper<Real, N, M>::DefaultH, typename ExplicitWrapper< Real, N, M >::FunctionJH Jh_x = ExplicitWrapper<Real, N, M>::DefaultJH, typename ExplicitWrapper< Real, N, M >::FunctionID in_domain = ExplicitWrapper<Real, N, M>::DefaultID )

inline

Set the explicit ODE/DAE system with lambda functions.

Parameters

[in]	name	The name of the explicit ODE/DAE system.
[in]	f	The explicit ODE system function.
[in]	Jf_x	The Jacobian of the explicit ODE system function with respect to the states.
[in]	h	The system's invariants.
[in]	Jh_x	The Jacobian of the system's invariants with respect to the states.
[in]	in_domain	The in-domain function.

explicit_system() [2/2]

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::explicit_system ( typename ExplicitWrapper< Real, N, M >::FunctionF f, typename ExplicitWrapper< Real, N, M >::FunctionJF Jf_x, typename ExplicitWrapper< Real, N, M >::FunctionH h = ExplicitWrapper<Real, N, M>::DefaultH, typename ExplicitWrapper< Real, N, M >::FunctionJH Jh_x = ExplicitWrapper<Real, N, M>::DefaultJH, typename ExplicitWrapper< Real, N, M >::FunctionID in_domain = ExplicitWrapper<Real, N, M>::DefaultID )

inline

Set the explicit ODE/DAE system with lambda functions.

Parameters

[in]	f	The explicit ODE system function.
[in]	Jf_x	The Jacobian of the explicit ODE system function with respect to the states.
[in]	h	The system's invariants.
[in]	Jh_x	The Jacobian of the system's invariants with respect to the states.
[in]	in_domain	The in-domain function.

has_system()

template<typename Real, Integer S, Integer N, Integer M>

bool Sandals::RungeKutta< Real, 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.

implicit_system() [1/2]

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::implicit_system ( std::string name, typename ImplicitWrapper< Real, N, M >::FunctionF F, typename ImplicitWrapper< Real, N, M >::FunctionJF JF_x, typename ImplicitWrapper< Real, N, M >::FunctionJF JF_x_dot, typename ImplicitWrapper< Real, N, M >::FunctionH h = ImplicitWrapper<Real, N, M>::DefaultH, typename ImplicitWrapper< Real, N, M >::FunctionJH Jh_x = ImplicitWrapper<Real, N, M>::DefaultJH, typename ImplicitWrapper< Real, N, M >::FunctionID in_domain = ImplicitWrapper<Real, N, M>::DefaultID )

inline

Set the implicit ODE/DAE system with lambda functions.

Parameters

[in]	name	The name of the implicit ODE/DAE system.
[in]	F	The ODE/DAE system function.
[in]	JF_x	The Jacobian of the ODE/DAE system function with respect to the states.
[in]	JF_x_dot	The Jacobian of the ODE/DAE system function with respect to the states derivative.
[in]	h	The system's invariants.
[in]	Jh_x	The Jacobian of the system's invariants with respect to the states.
[in]	in_domain	The in-domain function.

implicit_system() [2/2]

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::implicit_system ( typename ImplicitWrapper< Real, N, M >::FunctionF F, typename ImplicitWrapper< Real, N, M >::FunctionJF JF_x, typename ImplicitWrapper< Real, N, M >::FunctionJF JF_x_dot, typename ImplicitWrapper< Real, N, M >::FunctionH h = ImplicitWrapper<Real, N, M>::DefaultH, typename ImplicitWrapper< Real, N, M >::FunctionJH Jh_x = ImplicitWrapper<Real, N, M>::DefaultJH, typename ImplicitWrapper< Real, N, M >::FunctionID in_domain = ImplicitWrapper<Real, N, M>::DefaultID )

inline

Set the implicit ODE/DAE system with lambda functions.

Parameters

[in]	F	The ODE/DAE system function.
[in]	JF_x	The Jacobian of the ODE/DAE system function with respect to the states.
[in]	JF_x_dot	The Jacobian of the ODE/DAE system function with respect to the states derivative.
[in]	h	The system's invariants.
[in]	Jh_x	The Jacobian of the system's invariants with respect to the states.
[in]	in_domain	The in-domain function.

info() [1/2]

template<typename Real, Integer S, Integer N, Integer M>

std::string Sandals::RungeKutta< Real, 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<typename Real, Integer S, Integer N, Integer M>

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

inline

Print the Runge-Kutta method information on a stream.

Parameters

[in,out]

os

Output stream.

irk_function()

template<typename Real, Integer S, Integer N, Integer M>

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

inline

Compute the Jacobian of the variables \( \mathbf{K} \) with respect to the states \( \mathbf{x} \) as

\[ \frac{\partial\mathbf{K}}{\partial\mathbf{x}} = \begin{bmatrix} \mathbf{I} & \mathbf{0} & \cdots & \mathbf{0} \\ & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{I} & \cdots & \mathbf{0} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{I} \end{bmatrix} \text{,} \]

where \( \mathbf{I} \) is the identity matrix of size \( N \times N \).

Parameters

[out]

jac

The Jacobian of the variables \( \mathbf{K} \) with respect to the states \( \mathbf{x} \). 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]	x	States \( \mathbf{x} \).
[in]	t	Independent variable (or time) \( t \).
[in]	h	Advancing step \( h \).
[in]	K	Variables \( \tilde{\mathbf{K}} = h \mathbf{K} \) of the system to be solved.
[out]	fun	The residual of system to be solved.

irk_jacobian()

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::irk_jacobian ( VectorN const & x, Real const t, Real const 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]	x	States \( \mathbf{x} \).
[in]	t	Independent variable (or time) \( t \).
[in]	h	Advancing step \( h \).
[in]	K	Variables \( h \mathbf{K} \) of the system to be solved.
[out]	jac	The Jacobian of system to be solved.

irk_step()

template<typename Real, Integer S, Integer N, Integer M>

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

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_old	States \( \mathbf{x}_k \) at the \( k \)-th step.
[in]	t_old	Independent variable (or time) \( t_k \) at the \( k \)-th step.
[in]	h_old	Advancing step \( h_k \) at the \( k \)-th step.
[out]	x_new	Computed states \( \mathbf{x}_{k+1} \) at the \( (k+1) \)-th step.
[out]	h_new	The suggested step \( h_{k+1}^\star \) for the next advancing step.
[out]	K	The \( \mathbf{K} \) variables of the Runge-Kutta method.

Returns

True if the step is successfully computed, false otherwise.

is_dirk()

template<typename Real, Integer S, Integer N, Integer M>

bool Sandals::RungeKutta< Real, 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<typename Real, Integer S, Integer N, Integer M>

bool Sandals::RungeKutta< Real, 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<typename Real, Integer S, Integer N, Integer M>

bool Sandals::RungeKutta< Real, 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<typename Real, Integer S, Integer N, Integer M>

bool Sandals::RungeKutta< Real, 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.

linear_system() [1/2]

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::linear_system ( std::string name, typename LinearWrapper< Real, N, M >::FunctionE E, typename LinearWrapper< Real, N, M >::FunctionA A, typename LinearWrapper< Real, N, M >::FunctionB b, typename LinearWrapper< Real, N, M >::FunctionH h = LinearWrapper<Real, N, M>::DefaultH, typename LinearWrapper< Real, N, M >::FunctionJH Jh_x = LinearWrapper<Real, N, M>::DefaultJH, typename LinearWrapper< Real, N, M >::FunctionID in_domain = LinearWrapper<Real, N, M>::DefaultID )

inline

Set the linear ODE/DAE system with lambda functions.

Parameters

[in]	name	The name of the linear ODE/DAE system.
[in]	E	The mass matrix function.
[in]	A	The system matrix function.
[in]	b	The system vector function.
[in]	h	The invariants function.
[in]	Jh_x	The Jacobian of the invariants function.
[in]	in_domain	The in-domain function.

linear_system() [2/2]

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::linear_system ( typename LinearWrapper< Real, N, M >::FunctionE E, typename LinearWrapper< Real, N, M >::FunctionA A, typename LinearWrapper< Real, N, M >::FunctionB b, typename LinearWrapper< Real, N, M >::FunctionH h = LinearWrapper<Real, N, M>::DefaultH, typename LinearWrapper< Real, N, M >::FunctionJH Jh_x = LinearWrapper<Real, N, M>::DefaultJH, typename LinearWrapper< Real, N, M >::FunctionID in_domain = LinearWrapper<Real, N, M>::DefaultID )

inline

Set the linear ODE/DAE system with lambda functions.

Parameters

[in]	E	The mass matrix function.
[in]	A	The system matrix function.
[in]	b	The system vector function.
[in]	h	The invariants function.
[in]	Jh_x	The Jacobian of the invariants function.
[in]	in_domain	The in-domain function.

max_projection_iterations() [1/2]

template<typename Real, Integer S, Integer N, Integer M>

Integer & Sandals::RungeKutta< Real, 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<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::max_projection_iterations ( Integer const t_max_projection_iterations )

inline

Set the maximum number of projection iterations.

Parameters

[in]

t_max_projection_iterations

The maximum number of projection iterations.

max_safety_factor() [1/2]

template<typename Real, Integer S, Integer N, Integer M>

Real & Sandals::RungeKutta< Real, 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<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::max_safety_factor ( Real const t_max_safety_factor )

inline

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

Parameters

[in]

t_max_safety_factor

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

max_substeps() [1/2]

template<typename Real, Integer S, Integer N, Integer M>

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

inline

Get the maximum number of substeps.

Returns

The maximum number of substeps.

max_substeps() [2/2]

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::max_substeps ( Integer const t_max_substeps )

inline

Set the maximum number of substeps.

Parameters

[in]

t_max_substeps

The maximum number of substeps.

min_safety_factor() [1/2]

template<typename Real, Integer S, Integer N, Integer M>

Real & Sandals::RungeKutta< Real, 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<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::min_safety_factor ( Real const t_min_safety_factor )

inline

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

Parameters

[in]

t_min_safety_factor

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

min_step() [1/2]

template<typename Real, Integer S, Integer N, Integer M>

Real & Sandals::RungeKutta< Real, 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<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::min_step ( Real const t_min_step )

inline

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

Parameters

[in]

t_min_step

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

name()

template<typename Real, Integer S, Integer N, Integer M>

std::string Sandals::RungeKutta< Real, S, N, M >::name ( ) const

inline

Get the name of the Runge-Kutta method.

Returns

The name of the Runge-Kutta method.

operator=()

template<typename Real, Integer S, Integer N, Integer M>

RungeKutta & Sandals::RungeKutta< Real, S, N, M >::operator= ( RungeKutta< Real, S, N, M > const & )

delete

Assignment operator for the timer.

order()

template<typename Real, Integer S, Integer N, Integer M>

Integer Sandals::RungeKutta< Real, 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<typename Real, Integer S, Integer N, Integer M>

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

inline

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

Parameters

[in]	x	The initial guess for the states \(\widetilde{\mathbf{x}} \).
[in]	t	The independent variable (or time) \( t \) at which the states are evaluated.
[out]	x_projected	The 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<typename Real, Integer S, Integer N, Integer M>

bool Sandals::RungeKutta< Real, S, N, M >::project_ics ( VectorN const & x, Real const 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]	x	The initial guess for the states \(\widetilde{\mathbf{x}} \).
[in]	t	The independent variable (or time) \( t \) at which the states are evaluated.
[in]	projected_equations	The indices of the states to be projected.
[in]	projected_invariants	The indices of the invariants to be projected.
[out]	x_projected	The 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<typename Real, Integer S, Integer N, Integer M>

bool Sandals::RungeKutta< Real, S, N, M >::projection ( )

inline

Get projection mode.

Returns

The projection mode.

projection() [2/2]

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::projection ( bool t_projection )

inline

Set the projection mode.

Parameters

[in]

t_projection

The projection mode.

projection_tolerance() [1/2]

template<typename Real, Integer S, Integer N, Integer M>

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

inline

Get the projection tolerance.

Returns

The projection tolerance.

projection_tolerance() [2/2]

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::projection_tolerance ( Real const t_projection_tolerance )

inline

Set the projection tolerance.

Parameters

[in]

t_projection_tolerance

The projection tolerance.

relative_tolerance() [1/2]

template<typename Real, Integer S, Integer N, Integer M>

Real Sandals::RungeKutta< Real, 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<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::relative_tolerance ( Real const t_relative_tolerance )

inline

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

Parameters

[in]

t_relative_tolerance

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

reverse()

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::reverse ( bool t_reverse )

inline

Set the time reverse mode.

Parameters

[in]

t_reverse

The time reverse mode.

reverse_mode()

template<typename Real, Integer S, Integer N, Integer M>

bool Sandals::RungeKutta< Real, S, N, M >::reverse_mode ( )

inline

Get the time reverse mode.

Returns

The time reverse mode.

safety_factor() [1/2]

template<typename Real, Integer S, Integer N, Integer M>

Real & Sandals::RungeKutta< Real, 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<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::safety_factor ( Real const t_safety_factor )

inline

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

Parameters

[in]

t_safety_factor

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

semi_explicit_system() [1/2]

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::semi_explicit_system ( std::string name, typename SemiExplicitWrapper< Real, N, M >::FunctionA A, typename SemiExplicitWrapper< Real, N, M >::FunctionTA TA_x, typename SemiExplicitWrapper< Real, N, M >::FunctionB b, typename SemiExplicitWrapper< Real, N, M >::FunctionJB Jb_x, typename SemiExplicitWrapper< Real, N, M >::FunctionH h = SemiExplicitWrapper<Real, N, M>::DefaultH, typename SemiExplicitWrapper< Real, N, M >::FunctionJH Jh_x = SemiExplicitWrapper<Real, N, M>::DefaultJH, typename SemiExplicitWrapper< Real, N, M >::FunctionID in_domain = SemiExplicitWrapper<Real, N, M>::DefaultID )

inline

Set the semi-explicit ODE/DAE system with lambda functions.

Parameters

[in]	name	The name of the semi-explicit ODE/DAE system.
[in]	A	The function for the mass matrix.
[in]	TA_x	The function for the mass matrix.
[in]	b	The function for the right-hand-side.
[in]	Jb_x	The function for the right-hand-side Jacobian.
[in]	h	The invariants function.
[in]	Jh_x	The Jacobian of the invariants function.
[in]	in_domain	The in-domain function.

semi_explicit_system() [2/2]

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::semi_explicit_system ( typename SemiExplicitWrapper< Real, N, M >::FunctionA A, typename SemiExplicitWrapper< Real, N, M >::FunctionTA TA_x, typename SemiExplicitWrapper< Real, N, M >::FunctionB b, typename SemiExplicitWrapper< Real, N, M >::FunctionJB Jb_x, typename SemiExplicitWrapper< Real, N, M >::FunctionH h = SemiExplicitWrapper<Real, N, M>::DefaultH, typename SemiExplicitWrapper< Real, N, M >::FunctionJH Jh_x = SemiExplicitWrapper<Real, N, M>::DefaultJH, typename SemiExplicitWrapper< Real, N, M >::FunctionID in_domain = SemiExplicitWrapper<Real, N, M>::DefaultID )

inline

Set the semi-explicit ODE/DAE system with lambda functions.

Parameters

[in]	A	The function for the mass matrix.
[in]	TA_x	The function for the mass matrix.
[in]	b	The function for the right-hand-side.
[in]	Jb_x	The function for the right-hand-side Jacobian.
[in]	h	The invariants function.
[in]	Jh_x	The Jacobian of the invariants function.
[in]	in_domain	The in-domain function.

solve()

template<typename Real, Integer S, Integer N, Integer M>

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

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_mesh	Independent variable (or time) mesh \( \mathbf{t} \).
[in]	ics	Initial conditions \( \mathbf{x}(t = 0) \).
[out]	sol	The solution of the system over the mesh of independent variable.

Returns

True if the system is successfully solved, false otherwise.

Warning

Do not use the solution for internal backtracking, as the step callback may directly modify the solution.

stages()

template<typename Real, Integer S, Integer N, Integer M>

Integer Sandals::RungeKutta< Real, 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<typename Real, Integer S, Integer N, Integer M>

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

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_old	States \( \mathbf{x}_k \) at the \( k \)-th step.
[in]	t_old	Independent variable (or time) \( t_k \) at the \( k \)-th step.
[in]	h_old	Advancing step \( h_k \) at the \( k \)-th step.
[out]	x_new	Computed states \( \mathbf{x}_{k+1} \) at the \( (k+1) \)-th step.
[out]	h_new	The suggested step \( h_{k+1}^\star \) for the next advancing step.
[out]	K	The \( \mathbf{K} \) variables of the Runge-Kutta method.

Returns

True if the step is successfully computed, false otherwise.

step_callback() [1/2]

template<typename Real, Integer S, Integer N, Integer M>

FunctionSC Sandals::RungeKutta< Real, S, N, M >::step_callback ( )

inline

Get the step callback function.

Returns

The step callback function.

step_callback() [2/2]

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::step_callback ( FunctionSC const & t_step_callback )

inline

Set the step callback function.

Parameters

[in]

t_step_callback

The step callback function.

system() [1/2]

template<typename Real, Integer S, Integer N, Integer M>

System Sandals::RungeKutta< Real, S, N, M >::system ( )

inline

Get the ODE/DAE system pointer.

Returns

The ODE/DAE system pointer.

system() [2/2]

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::system ( System t_system )

inline

Set the ODE/DAE system pointer.

Parameters

[in]

t_system

The ODE/DAE system pointer.

tableau() [1/2]

template<typename Real, Integer S, Integer N, Integer M>

Tableau< Real, S > & Sandals::RungeKutta< Real, S, N, M >::tableau ( )

inline

Get the Butcher Tableau reference.

Returns

The Tableau reference.

tableau() [2/2]

template<typename Real, Integer S, Integer N, Integer M>

Tableau< Real, S > const & Sandals::RungeKutta< Real, S, N, M >::tableau ( ) const

inline

Get the Butcher Tableau const reference.

Returns

The Tableau const reference.

type()

template<typename Real, Integer S, Integer N, Integer M>

Type Sandals::RungeKutta< Real, 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<typename Real, Integer S, Integer N, Integer M>

bool Sandals::RungeKutta< Real, S, N, M >::verbose_mode ( )

inline

Get the verbose mode.

Returns

The verbose mode.

verbose_mode() [2/2]

template<typename Real, Integer S, Integer N, Integer M>

void Sandals::RungeKutta< Real, S, N, M >::verbose_mode ( bool t_verbose )

inline

Set the verbose mode.

Parameters

[in]

t_verbose

The verbose mode.

Member Data Documentation

m_absolute_tolerance

template<typename Real, Integer S, Integer N, Integer M>

Real Sandals::RungeKutta< Real, S, N, M >::m_absolute_tolerance {1e-6}

private

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

m_adaptive

template<typename Real, Integer S, Integer N, Integer M>

bool Sandals::RungeKutta< Real, S, N, M >::m_adaptive {true}

private

Adaptive step mode boolean.

m_lu

template<typename Real, Integer S, Integer N, Integer M>

Eigen::FullPivLU<MatrixP> Sandals::RungeKutta< Real, S, N, M >::m_lu

mutableprivate

LU decomposition for the projection matrix.

m_max_projection_iterations

template<typename Real, Integer S, Integer N, Integer M>

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

private

Maximum number of projection steps.

m_max_safety_factor

template<typename Real, Integer S, Integer N, Integer M>

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

private

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

m_max_substeps

template<typename Real, Integer S, Integer N, Integer M>

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

private

Maximum number of substeps.

m_min_safety_factor

template<typename Real, Integer S, Integer N, Integer M>

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

private

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

m_min_step

template<typename Real, Integer S, Integer N, Integer M>

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

private

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

m_newtonK

template<typename Real, Integer S, Integer N, Integer M>

NewtonK Sandals::RungeKutta< Real, S, N, M >::m_newtonK

mutableprivate

Newton solver for IRK methods.

m_newtonX

template<typename Real, Integer S, Integer N, Integer M>

NewtonX Sandals::RungeKutta< Real, S, N, M >::m_newtonX

mutableprivate

Newton solver for ERK and DIRK methods.

m_projection

template<typename Real, Integer S, Integer N, Integer M>

bool Sandals::RungeKutta< Real, S, N, M >::m_projection {true}

private

Adaptive step mode boolean.

m_projection_tolerance

template<typename Real, Integer S, Integer N, Integer M>

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

private

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

m_relative_tolerance

template<typename Real, Integer S, Integer N, Integer M>

Real Sandals::RungeKutta< Real, S, N, M >::m_relative_tolerance {1e-3}

private

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

m_reverse

template<typename Real, Integer S, Integer N, Integer M>

bool Sandals::RungeKutta< Real, S, N, M >::m_reverse {false}

private

Time reverse mode boolean.

m_safety_factor

template<typename Real, Integer S, Integer N, Integer M>

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

private

Safety factor for adaptive step \( f \).

m_step_callback

template<typename Real, Integer S, Integer N, Integer M>

FunctionSC Sandals::RungeKutta< Real, S, N, M >::m_step_callback {nullptr}

private

Step callback function.

m_system

template<typename Real, Integer S, Integer N, Integer M>

System Sandals::RungeKutta< Real, S, N, M >::m_system

private

ODE/DAE system object pointer.

m_tableau

template<typename Real, Integer S, Integer N, Integer M>

Tableau<Real, S> Sandals::RungeKutta< Real, S, N, M >::m_tableau

private

Butcher tableau of the Runge-Kutta method.

m_verbose

template<typename Real, Integer S, Integer N, Integer M>

bool Sandals::RungeKutta< Real, S, N, M >::m_verbose {false}

private

Verbose mode boolean.

SQRT_EPSILON

template<typename Real, Integer S, Integer N, Integer M>

const Real Sandals::RungeKutta< Real, S, N, M >::SQRT_EPSILON {std::sqrt(EPSILON)}

< Basic constants. Square root of machine epsilon epsilon static constant value.

---

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

• include/Sandals/RungeKutta.hh