COMMUNICATIONS IN NUMERICAL METHODS IN ENGINEERINGCommun. Numer. Meth. Engng (in press)Published online in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/cnm.963Dynamically equivalent implicit algorithms for the integrationof rigid body rotationsP. Krysl∗, †University of California, San Diego, 9500 Gilman Dr #0085, La Jolla, CA 92093-0085, U.S.A.SUMMARYTwo midpoint-trapezoid pairs of dynamically equivalent (conjugate) algorithms are derived as compositionsof first-order forward Euler and backward Euler integrators as applied to an incremental form of theinitial-value problem of three-dimensional rigid body rotation. The algorithms are related to the recentlydeveloped methodology of the so-called Munthe-Kaas Runge–Kutta methods. Selected examples are usedto illustrate the excellent long-term integration properties. Copyrightq2006 John Wiley & Sons, Ltd.Received 20 April 2006; Accepted 11 October 2006KEY WORDS: rigid body dynamics; time integration; trapezoidal rule; midpoint rule; symplectic;angular momentum conservation; dynamic equivalence; composition of maps; Munthe-Kaas Runge–Kutta methods1. INTRODUCTIONIt is well-known that for the time integration of vector-space Hamiltonian mechanics the midpointrule and the trapezoidal rule are dynamically equivalent [1]. The present paper addresses thequestion whether there is a corresponding pair of dynamically equivalent algorithms for the initial-value ordinary differential equation problem that describes three-dimensional rigid body rotations.The tool used in the construction of these algorithms is the composition of maps. It is commonknowledge that both the midpoint rule and the trapezoidal rule in the vector-space setting resultfrom the composition of the forward and backward (backward and forward, respectively) half-stepEuler methods. In the first part, this procedure is reviewed to provide a template for the secondpart of the manuscript. Following the template, not one, but two midpoint-trapezoid pairs ofdynamically equivalent (conjugate) algorithms are derived as compositions of first-order forwardEuler and backward Euler integrators as applied to an incremental form of the initial-value problem.∗Correspondence to: P. Krysl, University of California, San Diego, 9500 Gilman Dr #0085, La Jolla, CA 92093-0085,U.S.A.†E-mail: [email protected] John Wiley & Sons, Ltd.P. KRYSLThe properties of these algorithms are investigated using a few previously published integratorsas a reference. In particular, the algorithms are viewed in the context of the recently developedmethodology of the so-called Munthe-Kaas Runge–Kutta methods [2]. A few selected examplesare used to illustrate the excellent properties of the proposed algorithms, especially the stable andwell-behaved response for very long integration intervals.2. DYNAMICS ON VECTOR SPACESThe initial-value problem for a mechanical system (for instance, a system of interacting particles)described by a vector of configuration variables (displacements) u ∈ Rnmay be stated as˙p = f, p(0) = p0˙u = M−1p, u(0) = u0(1)where˙p is the rate of linear momentum,˙u is the velocity, and f = f(u, t) is the applied force. Forsimplicity we shall assume a time-independent mass matrix M. The initial values are p0,andu0.Using subscripts to indicate the time to which a given quantity belongs, and writing ft+h=f(ut+h, t + h), we may formulate a forward Euler time step applied to (1) aspt+h= pt+ hftut+h= ut+ hM−1pt(2)and the backward Euler approximation aspt+h= pt+ hft+hut+h= ut+ hM−1pt+h(3)Sequential composition of algorithms (2) and (3) in this order with time steps h = t/2 yields thetrapezoidal rule [1]pt+t= pt+t2(ft+ ft+t)ut+t= ut+t2M−1(pt+ pt+t)(4)Sequential composition of algorithms (3) and (2) with time steps h = t/2 yields the midpointrule [1]pt+t= pt+ tft+t/2ut+t= ut+t2M−1pt+t/2(5)The above algorithms have been derived as forward or backward Euler approximations to thederivatives in Equation (1), but, importantly, it would equally make sense to understand them asCopyright q 2006 John Wiley & Sons, Ltd. Commun. Numer. Meth. Engng (in press)DOI: 10.1002/cnmINTEGRATION OF RIGID BODY ROTATIONSapproximations of the integrals in these, equivalent, equationspt+h= pt+t+htf() dut+h= ut+ M−1t+htp() d(6)For the vector-space problem there is no advantage to be gained from either form, but as we shallsee next, these two approaches yield different algorithms for the rotational dynamics.3. DYNAMICS OF ROTATIONSAn excellent discussion of the difficulties of interpolating on the curved manifold SO(3), which isan appropriate setting for this problem, has been given in Reference [3]. To begin our presentation,we shall show how to formulate rotational dynamics algorithms to bring out the parallels to theabove vector-space trapezoidal/midpoint dynamically equivalent couple in Equations (4) and (5).The initial-value problem may be written in the convected description (body frame) as˙P =−skew[I−1P]P + T, P(0) = P0˙R = R skew[I−1P], R(0) = R0(7)where˙P is the rate of body-frame angular momentum, R is the rotation matrix (tensor), R−1= RT(orthogonal operator), T = T(t, R(t)) is the applied torque in the body frame, skew[•] is definedby skew[w]·w = 0,andI is the time-independent tensor of inertia in the body frame.3.1. Vector parametrizationThe second equation (7) is not in a form suitable for forward or backward Euler discretization:the rotation tensor constitutes points of the Lie group SO(3), which is not a vector space andlinear combinations are not legal operations on the rotation tensors. Therefore, an inevitable lossof orthogonality of the rotation tensor would result when the time stepping was applied directly.To transform the initial-value problem to a form suitable for our purposes, we shall introduce therotation vector representation o f the rotation tensor.As is standard, the equation of motion is written in the spatial frame as˙p = RT, where p = RPis the spatial angular momentum. Integrating the spatial equation of motion, and converting backto the body frame, we may write the equation of motion in integral form in the body frame asP(t + t ) = exp[−W]P(t) + RT(t)t+ttR()T() d (8)where exp[−W]= exp[−skew W]=RT(t +t)R(t) is the incremental rotation through vector (−W).By time differentiation of Equation (8), and by comparison with the original