|
| 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 |
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. |