Solving stiff multibody systems by implicit Runge-Kutta methods

Meike Schaub, Bernd Simeon

Centre for Scientific Computing and Mathematical Modelling (IWRMM), Universität Karlsruhe, D-76128 Karlsruhe, Germany

schaub@iwrmm.math.uni-karlsruhe.de

http://www.mathematik.uni-karlsruhe.de/~iwrmm/Persons/Schaub/english.html

Contributed talk

Multibody systems with stiff springs or discretized elastic bodies lead quite often to stiff differential equations as explicit numerical methods are forced to use very small time steps due to high frequency perturbations. To avoid this phenomenon, implicit methods such as implicit Runge-Kutta-methods are used which damp out high artificial frequency modes, see e.g. the stiff beam in [1].

But even for implicit methods standard stepsize strategies may yield small time steps due to the order reduction phenomenon. To avoid this, $h$-scaling is a convenient remedy as proposed in [3]. Hereby the error estimators of stiff velocity components are multiplied by the stepsize $h$ to overcome the reduced order, see [5] for examples. However, stiff components have to be known in advance and thus an analysis to locate them is advantageous.

In the talk different error estimators in [1] are investigated and an algorithm for stiffness detection based on the iteration matrix $(I-\gamma_0 h J)^{-1}$ is derived, where $J$ is the jacobian and $\gamma_0$ a constant. Implemented in combination with the codes RADAU5 of Hairer and Wanner [1] and SPARK3 of Jay [2] this approach shows a significant improvement of the stepsize sequence concerning smoothness and number of rejected steps. Details of the algorithm can be found in [4].

Another important issue of stiff mechanical systems is the control of numerical damping. For this purpose we present the class of blended Lobatto methods which is based on a convex combination of two methods in the Lobatto family, see [2]. By introducing a parameter $\theta$, additional flexibility is achieved to adapt the numerical damping behaviour to the needs of the system and single components in particular. For example a combination of Lobatto IIIA and IIIC leads to a class of L-stable methods that can be used to control the energy error in the linear case, the combination of the symplectic method Lobatto IIIAB and Lobatto IIIC yields a B-stable class of methods.