Explicit solvers

Explicit solvers are in general only applicable for systems without constraints (i.e., no joints!). However, some solvers accept simple CoordinateConstraint, e.g., fixing coordinates to the ground. Nevertheless, for constraint-free systems, e.g., with penalty constraints, can be solved for very high order and with great efficiency. A list of explicit solvers is available, see Section DynamicSolverType, for an overview of all implicit and explicit solvers.

The solution vector \(\txi\) (denoted as \(y\) in the literature ), which is defined as

\[\txi = [{\mathbf{q}}\tp \;\; \dot {\mathbf{q}}\tp \;\; {\mathbf{y}}\tp ]\tp\]

and which includes ODE2 coordinates and velocities and ODE1 coordinates. All coordinates are computed without reference values.

The ODE1 and ODE2 equations of Eq. (34), with \(\tlambda=0\), are written in explicit form and converted to first order equations,

(34)\[\begin{split}\dot {\mathbf{q}} &=& \vel \nonumber \\ \dot \vel & = &{\mathbf{M}}^{-1} {\mathbf{f}}_\SO({\mathbf{q}}, \vel, t) \nonumber \\ \dot {\mathbf{y}} & = &{\mathbf{f}}_\FO({\mathbf{y}}, t) \\\end{split}\]

The system first order differential equations for explicit solvers thus read

\[\dot \txi = {\mathbf{f}}_e (\txi, t)\]

Explicit Runge-Kutta method

Explicit time integration methods seek the solution \(\txi_{t+h}\) at time \(t+h\) for given initial value \(\txi_{t}\) (at the beginning of one step \(t\) or at the beginning of the simulation, \(t=0\)),

\[\txi_{t+h} = \txi_{t} + \Delta \txi.\]

For any given Runge-Kutta method, the integration of one step with step size \(h\) is performed by an approximation

(35)\[\Delta \txi = \int _{t}^{t+h}{\mathbf{f}}_e(\tau ,\txi(\tau ))d\tau \approx h\left[b_{1} {\mathbf{f}}_e(t,\txi(t))+b_{2} {\mathbf{f}}_e(t+c_{2} h,\txi(t+c_{2} h))+ \ldots +b_{s} {\mathbf{f}}_e(t+\txi_{s} h,u(t+\txi_{s} h))\right]\]

in which \(t + c_{i}h\) is the time for stage \(i\) and \(b_i\) the according weight given in the integration formula. Stages are within one step (therefor called one-step-methods), where \(c_i=0\) represents the beginning of the step and \(c_i=1\) the end. Note that \(c_{1}= 0\) for explicit integration formulas.

The unknown solution vectors \(\txi\) at the stages are abbreviated by

\[{\mathbf{g}}_{i} \approx \txi(t+c_{i} h)\]

and computed by explicit integration (quadrature) formulas of lower order (\(g_i\) not to be mixed up with algebraic equations!),

(36)\[\begin{split}\begin{array}{l} {{\mathbf{g}}_{1} =\txi_t} \\ {{\mathbf{g}}_{2} =\txi_t+ha_{21} {\mathbf{f}}_e(t,{\mathbf{g}}_{1} )} \\ {{\mathbf{g}}_{3} =\txi_t+h\left[a_{31} {\mathbf{f}}_e(t,{\mathbf{g}}_{1} )+a_{32} {\mathbf{f}}_e(t+c_{2} h,{\mathbf{g}}_{2} )\right]} \\ {{\rm \; \; \; \; \; \; }\vdots } \\ {{\mathbf{g}}_{s} =\txi_t+h\left[a_{s1} {\mathbf{f}}_e(t,{\mathbf{g}}_{1} )+a_{s2} {\mathbf{f}}_e(t+c_{2} h,{\mathbf{g}}_{2} )+ \ldots +a_{s,s-1} {\mathbf{f}}_e(t+c_{s-1} h,{\mathbf{g}}_{s-1} )\right]} \end{array}\end{split}\]

After all vectors \({\mathbf{g}}_i\) have been consecutively evaluated, the step is updated by Eq. (35).

For some exemplary tableaus of explicit and impliciti Runge-Kutta methods, see theDoc.pdf!

Automatic step size control

Advanced solvers, such as ODE23 and DOPRI5, include automatic step size control(activated with timeIntegration.automaticStepSize = True in simulationSettings).

We estimate the error of a time step with current step size \(h\) by using an embedded Runge-Kutta formula, which includes two approximations (35) of order \(p\) and \(\hat p = p-1\), which is obtained by using two different integration formulas with common coefficients \(c_i\), but two sets of weights \(b_i\) and \(\hat b_i\), leading to two approximations \(\txi\) and \(\hat \txi\). These so-called embedded Runge-Kutta formulas are widely used, for details see Hairer et al. .

The according apporximations \(\txi\) and \(\hat \txi\) are used to estimate an error

\[e_j=|\xi_j- \hat \xi_j|\]

for every component \(j\) of the solution vector \(\txi\). A scaling is used for every component of the solution vector, evaluating at the beginning (\(0\)) and end (\(1\)) of the time step:

\[s_j = a_{tol} + r_{tol} \cdot \mathrm{max}(|\xi_{0j}|, |\xi_{1j}|)\]

Then the relative, scaled, scalar error for the step, which needs to fulfill \(err \le 1\), is computed as

\[err = \sqrt{\frac 1 n \sum_{j=1}^n \left( \frac{\xi_{1j} - \hat \xi_{1j}}{s_j} \right)^2}\]

The optimal step size then reads

\[h_{opt} = h \cdot \left(\frac{1}{err} \right)^{(1/(q+1))}\]

Currently we use the suggested step size as

\[h_{new} = \mathrm{min}\left(h_{max}, \mathrm{min}\left(h \cdot f_{maxInc}, \mathrm{max}(h_{min}, f_{sfty} \cdot h_{opt}) \right) \right)\]

With the maximum step size \(h_{max} = \frac{t_{start} - t_{end}}{n_{steps}}\) and the minimum step size \(h_{min}\), given in the timeIntegrationsimulationSettings. The factor \(f_{maxInc}\) limits the increase of the current step size \(h\), the factor \(f_{sfty}\) is a safety factor for limiting the chosen step size relative to the optimal one in order to avoid frequent step rejections. If \(h_{new} \le h\), the current step is accepted, otherwise the step is recomputed with \(h_{new}\). For more details, see Hairer et al. .

Stability limit

Note that there are hard limitations for every explicit integration method regarding the step size. Especially for stiff systems (basically with high stiffness parameters and small masses, but also with restrictions to damping), the step size \(h\) has an upper limit: \(h < h_{lim}\). Above that limit the method is inherently unstable, which needs to be considered both for constant and automatic step size selection.

Explicit Lie group integrators

All explicit solvers including the automatic step size solvers (DOPRI5, ODE23) have been equiped with Lie group integration functionality, see Holzinger et al. .

Basically, the integration formulas, see Section Explicit Runge-Kutta method are extended for special rotation parameters. Lie group integration is currently only available for NodeRigidBodyRotVecLG used in ObjectRigidBody (3D rigid body). FFRFreducedOrder will be extended to such nodes in the near future. To get Lie group integrators running with rigid body models, all 3D node types need to be set to NodeRigidBodyRotVecLG and set explicitIntegration.useLieGroupIntegration = True.

Constraints with explicit solvers

Explicit solvers generally do not solve for algebraic constraints, except for very simple CoordinateConstraint. All connectors having the additional type=Constraint, see the according object in Section ObjectConnectorSpringDamperff., are in general not solvable by explicit solvers. Currently, only CoordinateConstraint with one coordinate fixed to ground can be accounted for, if explicitIntegration.eliminateConstraints == True. However, this offers the great flexibility to compute finite elements (imported meshes or ANCF beams) to be (partially) fixed to ground. A CoordinateConstraint that fixes a coordinate with index \(j\) to ground leads to the simple algebraic ODE2 equation

\[g_j({\mathbf{q}}) = 0 \quad \Leftrightarrow \quad q_j = 0\]

which can be solved by the implemented explicit solvers by just setting \(q_j = 0\) previously to every computation and \(\dot q_j = 0\) after every RHS evaluation.

NOTE that, if explicitIntegration.eliminateConstraints == False, constraints are ignored by the explicit solver (and all algebraic variables are set to zero). This may be wanted (e.g. to investigate the free motion of bodies), but in general leads to wrong and meaningless solution.