AP04-Bittnar1.vp 1 Introduction It is well-known that the Verlet algorithm (explicit New- mark for a certain choice of its parameters) may be written as a composition of two first order algorithms, the symplectic Euler and its adjoint [2]. What is perhaps less known is that there is an interpretation of this composition in terms of an approximation of the midpoint rule, which is of course im- plicit in the evaluation of the forcing impulse. Our goal in this brief note is to point out that this mechanically inspired deri- vation yields the well-known second order explicit Verlet algorithm in the vector space setting, and also an extremely accurate integrator for rigid body rotations when the mid- point rule is interpreted in the Lie sense. 2 Vector space midpoint rule approximation Let us write the initial value problem for a mechanical sys- tem with configuration u � �n in a fairly general form as � , ( ) � , ( ) p f p p u M p u u � � � �� 0 0 0 1 0, where �p is the rate of linear momentum, and f f u� ( , )t is the applied force. For simplicity we shall assume p Mv� , where M is a time-independent mass matrix, and v is the velocity. The midpoint approximation of the second equation is u u M p ( ) ( )t t t t t t� � � � � � � � �� � �1 2 , where p p u( ) ( ( ), )t t t t t t� � � �� � �2 2 2 makes this for- mula implicit. To approximate the midpoint momentum, we may use the equation of motion in integral form, p p f( ) ( ) ( )t h t t t h � � � � � � �d , and numerically evaluate the impulse integral. This we pro- pose to treat by recourse to the concept of concentrated, discrete impulses delivered at pre-selected time instants. In particular, the impulse may be delivered at the end of the time step, in which case f 0( )� �d t t t� � � 2 . On the other hand, the impulse may be imposed at the be- ginning of the time step, and f f( ) ( )� �d t t t t t � � � � 2 . Therefore, we obtain two algorithms. The first one, �� t , may be recognized as the symplectic Euler method, and the second, ��t * , as its adjoint. �� � � � �t t t t t t t t t t t t t p u p u p f u M p � � � � � � � � � � � � � � � � �1 � � � � � � � � t (symplectic Euler) �� � � �� � t t t t t t t t t t t t t * p u p u p f u M � � � � � � � � � � � � � � � � � � � � � � � �1p t (symplectic Euler adjoint). These algorithms are symplectic [2], and momentum con- serving. Their accuracy is only linear in the time step, but their composition preserves both symplecticity and momen- tum conservation and yields a second order accurate algo- rithm. That algorithm may be recognized as the well-known Verlet (explicit Newmark with � �1 2) � � �� � �t t t� 2 2 * � (Verlet). 3 Midpoint rule approximation on the rotation group Now we shall discuss the dynamics of rigid bodies rotating about a fixed point. (More detail is available in Reference [3].) The initial value problem may be written in the convected description (body frame) as � �� , ( )� � � � �� � � ��skew I T1 00 � �R R I R R� ��skew 1 00� , ( ) , where �� is the rate of body frame angular momentum, R is the rotation tensor, T is the applied torque in the body frame, skew[�] is defined by skew[ ]h h 0� � , and I is the time-inde- © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 51 Acta Polytechnica Vol. 44 No. 5 – 6/2004 Explicit Time Integrators for Nonlinear Dynamics Derived from the Midpoint Rule P. Krysl We address the design of time integrators for mechanical systems that are explicit in the forcing evaluations. Our starting point is the midpoint rule, either in the classical form for the vector space setting, or in the Lie form for the rotation group. By introducing discrete, concentrated impulses we can approximate the forcing impressed upon the system over the time step, and thus arrive at first-order integrators. These can then be composed to yield a second order integrator with very desirable properties: symplecticity and momentum conservation. Keywords: Time integration, rigid body motion, midpoint rule, symplectic Euler, Verlet, Newmark, midpoint Lie algorithm. pendent tensor of inertia in the body frame. The second equation is not in a form suitable for midpoint discretization, because the rotation tensor constitutes points of the Lie group SO(3), which is not a vector space and linear combinations are not legal operations on the rotation tensors. To transform the initial value problem to a form suitable for our purposes, we shall introduce the rotation vector representation of the rota- tion tensor. The equation of motion is written in the spatial frame as �� � RT, where �� is the rate of the spatial angular momen- tum, and consequently the equation of motion may be written in integral form as � �� � �( ) exp [ ] ( ) ( ) ( ) ( )t t tT t t � � � � � � � � � ��skew d0 0 R R T� � � � where � �exp [ ] ( ) ( )� �skew � R RT t t0 is the incremental rota- tion through vector �. Upon time differentiation and identi- fication with the original differential equation of motion, we obtain � ( exp )[ ]� ��� � � �d skew 1 1I , where d skewexp ( )[ ]� � � is the differential of the exponential map. The initial value problem may therefore be rewritten as � , ( )� � � � �� � � �� � �� � ��skew I T1 00 , � ( exp ) ( )[ ]� � ��� �� � �d ,skew 1 1 0I 0 . The midpoint approximation applied to the second equa- tion yields (�t � 0) � � � t t t tt t � � � � �� � � � � � �� � � d skew exp [ ] 1 2 1 1 2I which may be simplified by noting d skew exp [ ]� � � � �1 2 � � � t t t t t t � � � to give 1 2 � � � t t t t tI � �� �� . This equation needs to be solved for the rotation vector. Therefore, as for the vector space setting we get two differ- ent algorithms, depending on the chosen approximation of � t t� � 2. For the impulse applied at the beginning of the time step we obtain the SO(3) counterpart of the sym- plectic Euler integrator: � � � � � � � � t t t t t t t t t R R � � � � � � � � � � � � �� � �exp[ )]skew( ( ) exp[ )] � � t t t t t t�� � � � � � � � � T R skew( , where �t t�� solves 1 � �� � t tt t t t t tI T� � �� �� � � � � � � �� � �� �exp ( )skew 1 2 . On the other hand, the total torque impulse applied at the end of the time step yields the adjoint method � � � � � � � � t t t t t t t t t* exp[ ) R R � � � � � � � � � � � � �� � �skew( ]( ) exp[ )] � � t t t t t t t�� � � � � � � � � � � T R skew( , where �t t�� solves 1 � � � t t t t t tI � � �� �� � � � � � � �� � �� exp skew 1 2 . Both algorithms are first-order, symplectic, and momen- tum conserving. As before, the composition of these two algo- rithms in one time step provides us with a second order accurate algorithm, which is an analogy of the Verlet (ex- plicit Newmark with � �1 2) algorithm for the vector space dynamics � � �� � �t t t� 2 2 * � . 52 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No. 5 – 6/2004 Fig. 1: Fast Lagrangian top; on the left hand side convergence in the norm of the error in body-frame angular momentum; on the right hand side convergence in the norm of the error in the attitude matrix: AKW: implicit midpoint rule of Austin, Krishnaprasad, [1]: SW: Simo and Wong [5]; NMB: Krysl, Endres Newmark algorithm [4]; LIEMID[1]: implicit midpoint Lie: LIEMID[E1]: adjoint of the symplectic Euler (explicit midpoint Lie variant 1); LIEMID[E2]: symplectic Euler (explicit midpoint Lie variant 2); LIEMID[EA]: alternating explicit midpoint Lie method. The above algorithms have been called the explicit mid- point Lie variant 2 and 1 respectively, and the alternating midpoint Lie algorithm in Reference [3]. It bears emphasis that these algorithms are not simply the symplectic Euler and its adjoint. They all reduce to the full midpoint Lie algorithm for torque-free motion, and it is only the approximation of the midpoint momentum evaluation that distinguishes them from the implicit midpoint Lie rule. 4 Example The accuracy of the explicit midpoint Lie algorithms is rather remarkable, as may be seen in Fig. 1. We show conver- gence graphs for the fast spinning heavy top (the kinetic energy is dominant, and a numerical method has to deal effectively with precession and nutation which are motions of distinct frequencies). Even the first-order methods perform very well for larger time steps, and the present alternating midpoint Lie algorithm is the strongest performer out of a selection of the best currently available implicit and explicit algorithms, including the implicit midpoint Lie method. 5 Conclusions We have presented an approach to the construction of rigid body integrators, in particular for general 3-D rotations, that are explicit in the evaluation of the forcing, momen- tum-conserving, and symplectic. The starting point is the midpoint (implicit) rule, which is then treated with numerical integration of the forcing with concentrated impulses. The resulting algorithms conserve momentum, are symplectic, first-order, and they are adjoint. Consequently, their composi- tion leads to a second order algorithm, which may be readily interpreted as a Verlet (explicit Newmark) integrator, both in the vector space setting, and in the setting of the special or- thogonal group of rotations. (We would like to suggest as an appropriate name explicit midpoint Lie algorithm for the latter.) For rotational dynamics, this integrator is a new addition to the lineup of current high performance algorithms, and in fact numerical evidence suggests that it is the best second- -order integrator to date. Applications abound: molecular dynamics, micro magnetics, rigid body dynamics, finite ele- ment dynamics of deformable solids. Acknowledgment Support for this research by a Hughes-Christensen award is gratefully acknowledged. References [1] Austin M., Krishnaprasad P. S., Wang L. S.: “Almost Lie-Poisson integrators for the rigid body.” Journal of Computational Physics, Vol. 107 (1993), p. 105–117. [2] Hairer E., Lubich C., Wanner G.: “Geometric Numerical Integration. Structure-Preserving Algorithms for Ordi- nary Differential Equations.” Springer Series in Comput. Mathematics, Springer-Verlag, Vol. 31 (2002). [3] Krysl P.: “Explicit Momentum-conserving Integrator for Dynamics of Rigid Bodies Approximating the Midpoint Lie Algorithm.” International Journal for Numerical Methods in Engineering, (2004), submitted. [4] Krysl P., Endres L.: “Explicit Newmark/Verlet algorithm for time integration of the rotational dynamics of rigid bodies.” International Journal for Numerical Methods in En- gineering, (2004), submitted. [5] Simo J. C., Wong K. K.: “Unconditionally stable algo- rithms for rigid body dynamics that exactly preserve energy and momentum.” International Journal for Numer- ical Methods in Engineering, Vol. 31 (1991), p. 19–52. Petr Krysl e-mail: pkrysl@ucsd.edu University of California San Diego La Jolla California 92093-0085, USA © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 53 Acta Polytechnica Vol. 44 No. 5 – 6/2004