Computer Physics Communications ELSEVIER
Computer Physics Communications 98 (1996) 153166
Numerical simulation of a rigid rotating body by Obrechkoff integration David P. Stapleton Department of Mathematics and Statistics, University of Central Oklahoma, Edmond, OK 73034, USA
Received 22 January 1996
Abstract It is first supposed that Euler's differential equations for the rotational dynamics of a rigid body have been solved, so that the angular velocity vector for the body and the derivative thereof are known vectorvalued functions of time in a bodyfixed coordinate system. A method is discussed by which a set of ordinary differential equations (ODEs) may be numerically integrated to deduce the orientation of the body in an inertial frame. The method employs a quaternion formulation for the orientation of the body, exhibits a local truncation error which is typically of fifth order, adheres to the boundary of the region of absolute stability, and requires only one angular velocity evaluation and one angular acceleration evaluation each step. It is then supposed that the solution to Euler's equations for a given body is not known. An iterative process is obtained for the integration of Euler's equations and the ODEs for body orientation simultaneously. Keywords: Rotational dynamics; Angular momentum; Euler's equations; Rigid body; Satellite dynamics; Quatemions
1. A review of the quaternion formulation for rigid body dynamics Given a rigid rotating body, consider a coordinate transformation between a body frame of reference in which coordinates x b, Yb, Zb correspond to axes fixed in the body and an inertial frame of reference in which coordinates x i, yi, zi correspond to nonrotating axes whose origin coincides with the origin o f the body axes. This can be accomplished with the matrix of direction cosines G, in which each component Cij = C u ( t ) is a function of time and denotes the cosine o f the angle between the ith body axis and the j t h inertial axis, as
Yu ~/C21
C22
c2311yi.
LC3, c32 c,,JL , In the quaternion development, socalled " E u l e r p a r a m e t e r s " a, b, c, d are introduced such that 4a 2 =
1 + Cll
4 b 2 = 1 + C22
 C22 
C33 ,
C33 
Cll ,

00104655/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved. PII S00104655(96)000707
(1)
154
D.P. Stapleton/ Computer Physics Communications 98 (1996) 153166 4 c 2 = 1 + C 3 3  Cll  C22, 4 d 2 = 1 + Cll + C22 + C 3 3 , 4 a d = C23  C32, 4 b d = C31  C13, 4 c d = Ct2  C2j,
and C is parameterized in the form a 2 + d 2  b 2  C2 C =
2(abcd)
2 ( a b + cd) b2+d2a2c
2( ac + bd)
2( bc  ad)
2( ac  bd) 2
2(bc+ad)
(2)
c2 + d 2  a2  b E
By use of the notation x = [ a b c d] T, the Euler parameters can be shown to satisfy
I xl
= 1/8 2 + b 2 + c 2 + d E = 1,
(3)
and also the system of differential equations, 2 =
I2x,
(4)
where the dot denotes differentiation with respect to time t and O is a 4 × 4 matrix given in terms of roll, pitch and yaw rates (components of angular velocity in the x b, Yb, Zb frame) p = p ( t ) , q = q(t), r = r ( t ) by
o r 1 ~Q=2
r q p
0 p q
p 0 r
(5) "
If p , q, r are given C s functions on an interval [t o, tf] then standard theorems from differential equations show that (4) together with arbitrary initial values a(0) = a o, b(0) = b 0, c(0) = c 0, d(0) = d o has a unique solution which is at least C 6 on [t 0, tf]. I f the body axes and inertial axes are to coincide at t = 0 then the initial values a(0) = 0, b ( 0 ) = 0, c ( 0 ) = 0, d ( 0 ) = 1 are applied. Hence, i f E u l e r ' s equations have been solved to obtain C 5 roll, pitch, and y a w rates p, q, r f o r a rotating rigid body then the orientation o f the body as a function o f time may be deduced in terms o f (1) and (2), f r o m the C 6 Euler parameters a, b, c, d obtained as the solution o f the ODEs (4) f o r given initial conditions ao, bo,c o, d o. Eq. (3) holds identically for any solution o f (4) whose initial condition satisfies (3) and insures that the matrix C in (2) is a rotation matrix; however, numerical solutions typically fail to satisfy (3).
2. The variable Obrechkoff method for numerical integration It is desired to obtain a numerical solution for a system of the form (4)   subject to some initial values a 0, b 0, c 0, d o that satisfy (3), such that (3) holds at each step. In this section a suitable method is introduced. For convenience of notation it is first written for the classical (single) ODE d y / d t = f ( t , y ) and then it is applied to the system (4).
i The development of this material (with a few typos) may be found in [1].
D.P. Stapleton / Computer Physics Communications 98 (1996) 153166
155
The following notation is used to numerically integrate d y / d t = f ( t , y) with y(t o) = Yo on [t 0, t f ] in K steps: y(t) is the analytical solution on [t 0, tf], t I < t 2 < t 3 < • • • < t K_ I are successive times in (t 0, tf), n = any element from {0, 1. . . . . K  1} "is the step number, h = t.+ ~ t. is the current step size. y. is the numerical estimate of y(t.), f . = f ( t . , y.), and f.~ is the total derivative.
af
af
f" = ~t(t,. y,) + ~y(t,. y,) f(td. y.). An implicit lstep Obrechkoff method 2 for d y / d t = f ( t , y) when f is differentiable is
y.+ , = y. + h( olf. +/3f.+ t) + h 2 ( Y f " + ~$f~+ l),
(6)
where a , /3, y. 6 are constants. This suggests the definition below, which may be extended in an obvious manner in order to integrate systems d y / d t = f ( t , y).
Definition. Let f be differentiable on a region D={(t,y)
lt o < t < t
e,  ~ < y < ~ c } ,
and let Yo be any real number. A numerical integration of the form (6) for d y / d t = f ( t , y), y(O) = Yo in which some or all of a = or(h), /3 = / 3 ( h ) , T = y(h), and 6 = 6(h) depend upon h shall be called an implicit 1step variable Obrechkoffmethod for d y / d t = f ( t , y) and the corresponding coefficients a , /3, Y, 6 shall be called
variable coefficients. If f is at least C 5 on D then substitution of the usual Maclaurin series expansions, Y n + l = Y n + h f n + 2 " J n±t.2r, + 6 " a n ±z.3¢,,+ . . . ,
I h2/w f.+ 1 = f . +hf" +7.. Jn
]
" "" ,
~ f,]+l =f~p +hf" z l n ' q '  2! l,~, J2 1h: m
"4"'',
into (6) shows that y . + , = y. + h( a + 13)f. + h2( 13 + y + 6 ) f " + h3(½/3 + 6 ) £ ' + h5('~4/3 + g v / j n
+ O(h6).
+
h4(//3
I ½8 ) / n "
(7)
Hence, for fifth order local truncation error it is required that u + 13 = 1 + v(h),/3 + y + 6 = ~ ' + g(h),½/3 + 6 = l + ~i = ~ + re(h), where v(h) = O(h4), g(h) = O(h3), k(h) = O(h2), re(h) = O(h). These gI + k ( h ) , and 3/3 equations determine that
or= 21 + v( h)  2 u ( h ) , 1
(8a)
fl=2+2u(h),
(8b)
Y= Wz + g( h)  k( h)  u( h),
(8c)
=

+ k(h)

u(h),
2 For a short description of these methods see [2].
(8d)
156
D.P. Stapleton / Computer Physics Communications 98 (1996) 153166
where u(h) = 3k(h)  6m(h) = O(h). The stability polynomial associated with (6), for the usual test equation d y / d t = Ay (A = a complex constant), is p( r) = r  1  h( aA + rflA)  h2(yA 2 + r3A2).
(9)
If h is such that the root r of p ( r ) satisfies ] r l < 1 then the method is absolutely stable and if h is such that the root of p ( r ) satisfies I r ] = 1 then the method has obtained the boundary of the region of absolute stability (for the given ODE and step size). If (6) is applied to a linear system, N
Yi= ~ aii( t) Yj, j=l
it is required that for each eigenvalue A of the matrix [a U] the root of the corresponding polynomial p ( r ) must satisfy I r I < 1 for absolute stability, or ] r I < 1 with I r [ = 1 for some eigenvalue to obtain the boundary of the region of absolute stability. Now let the implicit 1step variable Obrechkoff method be applied to ~: = 12x, and define ~=
1 ~p,
1 q = ~q,
~=[p~]~
~= j1r ,
,~2=,~,~, i,~1=¢~=#2+~2+~2
Let x . = [ a. b n c. d.] T denote the numerical estimate of x ( t . ) = [a(t.) b( t.) c( t.) d( t.)] T and 12. = 12(t.). The eigenvalues A+ = _+i I &l (each of multiplicity two) of 12 yield corresponding roots. (1
h2to2"y) +ihl &[ ct
r + = (I + h2t~2t$) Tih] &l f l '
(10)
of (9) so that absolute stability occurs if [ r+ ] < 1 and the boundary of the region of absolute stability is obtained if ]r+ } = 1. If fifth order local truncation error is required then Eqs. (8) show that ] r+ I < 1 is equivalent to 2 g ( h ) + v ( h ) [ 4 u ( h )  v ( h )  1] + h2wZ[g(h)  2 u ( h ) ] [ 2 k ( h )  g ( h )  6] > 0.
(11)
One way to obtain equality in (11) is to choose v(h) = g(h) = u(h) = O. Thus, the boundary of the region of absolute stability and at least fifth order local truncation error is attained if the method (6) with variable coefficients is applied to the system of ODEs (4) with
~=t3= ~, v =  ~ = ~ + k ( h ) ,
(12)
1
where k(h) = O(h 2) is arbitrary. In matrix form this method is written implicitly from (6) and (4) as
[, ~h12.+, where 6 =
12
~h2( a.+, + 12:+,)]x.+, = [, + ~h12. ~h2( a. + a~. )]x..
k(h). It can be rewritten using the fact that 12 2 = _ &2i, the notation
~2 = (1 + h~2)2 + (~h~ + h28~)2 + (~h~ + h2~)2 + (~h~ + h2~)2 , and the matrix inverse (for ~2:96 0) 1
1
157
D.P. Stapleton / Computer Physics Communications 98 (1996) 153166
in the explicit form (in which it is assumed that g2.+ 1, 1"2.+ ~, /z. ^ ~+ 1 and ton+^ 2 [ are known)
x.+l=
[(1 +hZTdo~+t)l+½hg2.+l
^2
+hZ3g'/.+l][(1 + h Z 3 & Z ) l + ½ h 1 2 .  h Z 3 ~ 2 . ] x . .
/~n+ 1
I
(13)
According to (10) and (12), this integration achieves the boundary of the region of absolute stability for any number 6. Also it has been shown that if 3 is the variable coefficient 3 = 12 k(h) then the local truncation error is O(hS). In terms of components, (13) is an+ 1 bn+ I Cn+ I d.+ 1 1
1 + h2&5~+ 1 1
2
^
~hr.+
1 ^2 P'n+ 1
1
l  h 3r.+
^
2

1 ^ ihq~  h 2 ~q.A 1
^
2
 2hPn+
l
1
^
2hq.+
2hr: ^  h28~,
½brAn + hE~f',
~
1
½h~,
1 + h23~02
1
^
2
2
l  h 3 q . + 1
~"
 ~hq. + h 3 q n
^
2
1 ^  ~hr, + h 2 3 rA.
^
2
1
^
2
~hq,+ l ÷ h 3q.+ 1
1 q h23t.~2+ 1
1 ^ 1 + h23~. + 1 ~hrn+
I
1 + h23&~
1
2h~.+ 1 at h 6p.+ l
1 ^ ~' ~h~.+ 1 + h 23p.+ 1
^
 ~hrn+
 hE3~n
,~
~hq.+ 1  h 3qn + 1
1
 2hq: ^ + h23~,
  ~1h ~^ + h23~.
 x h ~ n + h 6 p .

X
 ,i h P^ . + 1 _ h23~.+
X
A
1 + h23Co2.+
l
2hqn+ 1 + h 3q.+ 1
1 + h23& 2
2
^
~hr.+ 1 + h 3r.+ l
1
h
2
^2 + 1 1 + h 23to.
A
8rn+
1
£h~: ^  hE3pnX ~hq,l ^ 1 ^ ~hr.
1 +
_ h234.
h23~n
an1 b,
(14)
Cn
d.
h28Co2
Function k ( h ) = O(h 2) is arbitrary in (12) which suggests that 6 might be chosen to attempt to Eq. (3) holds for arbitrary initial condition x 0 such that I x01 = 1 if and only if (14) is an transformation. Since the first matrix in (14) has orthogonal columns of length /2.+ I and the second orthogonal columns of equal length the transformation is orthogonal if and only if the columns of matrix are of length t2,+ 1 This requirement is tantamount to
tOn+ltOn
O.)n+lO~ n + 3
2 ¢ ~ 2 + l   2 ^ton 2+ht~n+l
A
" tx ° n + l I h ~ n
" ton
"] X
satisfy (3). orthogonal matrix has the second
tOn+ l   COn
=
(15) where
and is a quadratic equation
AB2 + B 3 + C = O ,
(16)
D.P. Stapleton / Computer Physics Communications 98 (1996) 153166
158
in which A, B, C depend upon h, ~ . , &n+ i, ~ . , ('On+ ~ 1" Some insight into (15) is gained from ^ x 1 3X h&. = mn+l  & n  71 h 2ton~h ton +
O(h 4)
A ^ ^ I r2x 1 3".A~ hoJn+ 1 = tOn+ I   ton I1 ~/'/ ('On÷ i  ~ h ¢0n t O ( h 4 ) o
These relations imply that h2((.On+,'~2 [ __ (.~92)= (h£~)n+l ,I h ~ o n ) . (h(~)n+ 1  h ~ ) n )
[ 2 ( t o^n + l   t O n )^" l   : ; f
=
^ = [2(tO.+,
I h 2 ( ¢"On x + 1 _~on
)
' 3"~: 3X + tOn , h,o.
_ _ 51h t3X o.+
O(h4)] . [lh2(~ tOn+l +to. R ) + O ( h 4) ]
O(h4)].[,h2(2, n
_~.h'~On)_l_O(h4)]
= h 2 ( ('°n+ l  £ ° n ) " ( 2 ~ ) n "1 h'~)n) + O ( h 5 ) ,
and ^ h ~n+l"
~ ~n+l
A ) + ~n" ~n
"~" ("On+ l  ton
2 h2
^
.)
x
 ~0n" ~0
1 . 3 [ t ton+ ^  "~n
l
°)
"" F t'O. " "ff) + O ( h " ['On
4)
X ton X + ~'On" "(o.)  31h t3o . . ~on + O(h 4) = m^2. + l  o g^2 . + T n It3/ ~ to." 1 3A "" l 3^ "" ^2  to. ^1 + 7ho.~ . ~,n + gh = W.+l w n" ~ . + O ( h 4 )
(terms not shown contain higher derivatives at t = t n) so that by substitution into (15),
o).+,
~.)"
+ 6 [ 3 ( & ~ + ,  m^ n2 ) +
h 3 l7~
2~.+h~.
x + .,'to,,
h 3 lgt ~
+ O ( h 3)
n ' &"'°n + O ( h 4 )
]
1 ^2 +'i(Wn+,&2) =0"
(17)
In the special case d(&2)/dt 4= 0 at tn (i.e. &," ~ , v~O), when h is small enough, (17) can be written in the
form (16) with =
ton+l
l A
On
R
^2 ^2 ('On + 1  ('On I ^
+ O(h4)'
"'°
B = 3 + h 3 "fit'On " ton + ~ton " COn + O ( h 3 ) , ^2 ^2 On+ 1  ton
C=i.
So, by t^o . +
1 1  ton =
A = h212 ^2 ,on
ht~.
+ O(h2),
~ 1 =~.+h~.+O(h2),and ton+
t~2n +
^2 = 1   On
(~°n'~°nl]+O(h3); B = 3+h 2(°"'f°n+3~'°n'~°n +O(h ton^JJ
2h& . ~ . + O ( h 2 ) ,
3)
'
C=~.
D.P. Stapleton/ ComputerPhysics Communications98 (1996) 153166
159
If h is small then B > 0 and B 2  4AC > 0, so the root of smaller magnitude of (16) can be written
(
)
1 I+¢I_4AC/B
8= (~~)
(18)
2
1 + O ( e 2) for (18) with e = 4 A C / B 2 ( = O(h2)) gives The Maclaurin formula 1/(1 + 1/l  e ) = ~l + ~e
2c(
ac) + O ( h
= mn
=
} 4 " ~
(
c AC2 n 3 +O(h
4)
n
)
4)
1
12+h 2 n
• to. + 3 to." &.
+ O ( h 3)
• ton
 ll+
h 2 (¢5.~J.+2~..~.3&. • ton
432
2,~+
~
tO to.
+O(h 3)
2,o.~) + O(h~).
This result corresponds to h 2 ( iS. ~ . + 2t~n • ~ n
k( h) =  ~
3&. . ~o.
2~) +O(h~) I
in (12) so that substitution into Eq. (7) with a = fl = 2 and 8 =  y = the smaller magnitude root of (15) identifies the local truncation error in (13) for the special case d(t~E)/dt ~ 0 at t. as
x*
n+l   X " + I
= h5 
2tO2
3&,, to. x
4~
1 _(5)]] + O ( h 6 ) • X."['7"~"~'n
(19)
•
Here x * ( t ) denotes the result of exact integration from t. to t.+ I, x~+ l = X * ( t . + l ) , and x~5) denotes d S x * / d t 5 at t.. Alternatively, suppose that at t., d ( ~ E ) / d t = 0 but d2(t~2)/dt 2 4= 0 (i.e. &.. t~. = 0 but ~ . . tS. + r~2 ~ 0). Again (17) may be divided by to. ^ 2+ l  to.^2 for small h, and the relations
=
h~n 'F ~ n
tO n
"st~n tOn + O ( h 4
x
I h2X
,^ h 2 ( xt o .  t ^o . + t ox2) . + h 3(xt o n . txo . + ~ w .  t o ."x) + O ( h 4 ) ,
x2  ton x2 ~ 2 h ~ . tOn+l
t~.+h2~2+O(h3),
imply that
A=h
2 (~n • ton ~ ^
ton
x
• ton
+
A2
ton
+ O ( h 2)
B=3+~h
^ n • tO "A*n I 3 to A n • r~ "" n tO t~ n
.
/'2
x n ].tO n tO
+ O ( h 3)
C =±
160
D.P. Stapleton/ ComputerPhysics Communications98 (1996) 153166
in (16). The root of smaller m a g n i t u d e can be rewritten from a Maclaurin e x p a n s i o n as
C 8=
35 ,~ x h Co,, &. + T~tO. " to.
AC 2
B
"'"
B 3 +O(h2)=~+
18
t~.to.X +to.,~2
+O(h2)'
and does not satisfy (12); so the resulting local truncation error (for this case) is only ^
"A"
35 ,~
O(h4)
X,
X* tOn " tOn d "i~tOn " tOn .+ t  X.+ 1 =  "~8h4 .~. k O ( h S ) . t3. • ~ . + t~O n2
(20)
Finally, in the case t5 2 = a constant on [t.,t.+ j], 8 = 0 is the smallest m a g n i t u d e root in (15) and does not satisfy (12). F r o m (7) it is observed that in general the error is only O(h3), i.e.
I x;+l  Xn+l = ~2h3xn + O ( h4). [
(21)
Thus, the following theorem has been established for the rotating rigid body problem (4).
Theorem. C o n s i d e r an initial value p r o b l e m ~ = O x , x ( t . ) = x . in which x = [ a b c d] T denotes a vectorvalued function of t, x . is a unit vector, and /3 = / 3 ( t ) , ~ = ~ ( t ) , P = P(t) are given C 5 functions on [ t . , t.+ t] which determine
0
/3] 0 /3
°=
/3 0
/3
o
&. = [/3(t.) ~(t.) P(t.)] T,
to.^Z&.. &.,
to.=
~2 ~ ton • ton
p ( t . ) q( t.) P( t.
,
ton.
Let h = t.+ l  t.. If either (i) (ii) (iii)
d ( & 2) d~ :~ 0
dt
0
at t. (i.e. &. • t~. ~ 0 ) , and
  4 = 0 dt 2
at t . ~ti . e , t^o . . t o . = O a n d
to.&.+to.
]
or
t5 2 = a constant on [ t. t. + l ],
then for h sufficiently small, (15) has real solutions and when the root 8 o f smallest m a g n i t u d e is applied in the numerical integration (13) (or (14)) the local truncation error is (19), (20), or (21), respectively 3. Also, [ Xn+ 1 [ = 1 and the boundary o f the region of absolute stability is attained.
3 Caution is advised if d(to2)/dt = 0 at tn but equality does not hold; the analysis for very small h (in the derivation of(19)) shows that the error is O(hS), however, close consideration of the development of (19) and (20) shows that for somewhat larger values of h the O(h 4) term in (20) can better approximate the actual error. Similar results hold for borderline cases of types (ii) and (iii).
D.P. Stapleton / Computer Physics Communications 98 (1996) 153166
161
The lack of a formula f o r the truncation error in cases other than (i), (ii), (iii) o f the theorem suggests the value 6 =  ~ f o r these cases. This does not precisely preserve a 2n + b n2 + c 2n + d n2 = 1 in general (although an exception occurs if & is a constant vector); however, the truncation error is
[ Xn+• t  Xn+ 1 ~ "7fO , hs o) + O(h6) I °* n
(22)
and, according to the comment following (13), the method clings to the boundary of the region of absolute stability. It is also suggested that h be chosen small enough that
max(hll&_
h&~l) << I1 + h 2 ~
2]
+
"2 on [t n, tn+ 1], to ensure an accurate calculation of/xn+ l 4= 0 and so that matrices in (14) are strongly diagonally dominant. Care shouM be taken in calculating the roots of (15) since the quadratic formula suffers from substantial rounding errors for B z >> 41 AC I. One difficulty with the theorem is that it is not clear how small h must be chosen in cases (i) and (ii) for (15) to have real roots and for error estimates to be valid. The following strategy is reasonable for choosing h and estimating the local truncation error in cases (i) and (ii):  from an initial estimate h of the step size, halve h until Eq. (15) has real roots with B > 0, and the matrices in (14) are strongly diagonally dominant,  let 6 be (the root of smaller magnitude of (15)) computed by (18),  integrate with (14), and  estimate the error by (7), using 6 =  ~2 + O(h2) or 6 =  ~2 + O(h), so that respectively
.   X n + l = "7~'~ l "~'n Xn+¿
h3(
+
+ O(h6),
or
x:+ l  xn+ = ha(T2 + ~))?n + O ( h 5 ) .
3. Integration when the solution to Euler's equations is not known In a body coordinate system x b, Yb, Zb that uses principal axes, Euler's equations are IlP
( I 2  I 3 ) q r = ~'l,
1 4 (
6)rp =
13 r  ( Ii  I 2) Pq = r3, where I~, 12, 13 are the moments of inertia and %, ~2, ~3 are the components of the moment about the Xb, Yb, Zb axes. These can be written in terms of J1=2(1213)/11,
L=rJ211,
J2=2(13I1)/12,
M=~2/212,
N=r3/213,
J3=2(1112)/13,
162
D.P. Stapleton / Computer Physics Communications 98 (1996) 153166
and p, 4, f as
p =Jl4~+L, ~ = J2 ~/3 + M ,
(23)
= J3/34 + N. The Obrechkoff technique (6) with variable parameters (12) applied to (23) gives
I
Pn+ I qn+ 1 ~n+ l
= q. +½h J2(r~.+P.+l~.+l)+M.+M.+l ~.
[J3(p.4.+~.+14.+l)+N.+N.+ , ^2  4 ~ . ) + J 2 rn+,p.+,r.p. Jl J3 q.+l/3.+l +4n+tNn+,OnNn +~n+lMn+,~nMn] "
__
+L.+,L. "2^
+h2~
+?.+,L.+I?.L.+p.+IN.+ 1p.N.] +M.+I  M . J3 J2( ^2pn+lrn+ ^    P n r . ) + J l 4~+1rn+142~.
+~.+ i g . + ,  ~ . g . + 4.+,L.+, 4.L.] +U.+,U. which can also be written as ^
Pn+ l =
+ ~[,,(~.~. + ~.+,.~.+,)+'. +'.+,1 ^
qn+
l =
[1 h2~J2(J,~;+, + J,p.~+,)]'{4.[1  h2~J2(J,~.~ + J,~.~)] +½h[ J2(?./3. + r~.+, ~.+,) + M. + M.+,]
^
gn+ l
( P.q.+P.÷lq.+l ^ ^ ) +N.+N.+I +½h[ J3
]
(24)
D.P. Stapleton/ ComputerPhysicsCommunications98 (1996)153166
163
These equations together with (14) may be used iteratively each step to try to converge to a variable Obrechkoff solution of the form [Pn+ 1 qn+ l rn+l a~+~ b,+ l c,+ i d,+ l]T of Euler's equations (23) and Eq. (4) for orientation simultaneously. A prediction of [/3,+ l qn+ t ?n+ l] T, e.g.
1 inl 'lnn+L1
:~ : = ~n + hl Jz?nPn + M. L~n+IJ
LrnJ
(25)
+ ½h2l Jz(Jlr^2qn + j3p2qn + ?nLn + p.Nn) + l(In
LJ3.O.+un
is needed to estimate quantities on the right of (24) for the first iteration 4.
4. Examples
To illustrate the variable Obrechkoff method, consider integration for the orientation of a body which maintains constant, to = [p q r] T = [4 2f3 6]T with x 0 = [a 0 b 0 c o d0]T = [1 0 0 0]T using constant step size h = g t (a constant step size is applied only for comparison purposes and generally a variable step size is preferred). Orientation may be deduced by the rotation matrix (3) from the analytical solution of (4), [~1 a
1[
4 c°s(nt) ]  3 sin(4t) 7~'sin(at)  2 sin(4t)
The result a = cos(4t) is compared with the Obrechkoff estimate (14) using 6 =  T2 (since to is a constant 1 vector) and the classical fourstage RungeKutta estimate in Fig. 1, for h = 7. With this constant step size each numerical method has fourth order global truncation error, i.e. fifth order local truncation error, so that accuracies at early times (e.g. 010 sec.) are comparable. As t increases (e.g. 90100 sec.), however, the Obrechkoff estimate satisfies (3) and maintains the correct amplitude while the RungeKutta estimate fails to satisfy (3) and decays. The superior stability characteristics of the Obrechkoff integration appear to cause a slower divergence from the true solution (see the plot for 390400 sec.) also. Similar results occur for parameters b, c, d. The Obrechkoff estimate also has a local error estimate obtained from (4), (22), and 0 2= ~21, of
x ; + ,   x . + t = 2~6hSO.xn+O(h 6) = 1.04)<
l0 6
[03 3
2
0 2
2 0
~ 3
bn
f3"
3
0
d,
c,
"
Next. consider o J = [ p q r ] v = [ 4 4 t 2 t ] r with x o = [ a o b o co do]T = [ 1 0 0 0 ] T and constant step size h = ~ . Fig. 2 compares the Obrechkoff estimate (14), using 6 computed by (18) when possible for h = ~ and
4 The predictorcorrectorschemeof (25) followedby only one correctionwith(24) is not suggestedsince(3) willnot hold in generaland because the stabilityproperties of such a methodare suspect.
~
II
8"
~j
~>
°
ii
ii
iir
i i Lir
.o ........
i ,,
~0 .......
i _ l l ~ L
~  ~ ~ _~ ~~
iii
iii
LLt
i Fi
i i t ~
...... "
bl
I~.~ ll
o
t
t
t
t
P~al~e~er
a
D.P. Stapleton / Computer Physics Communications 98 (1996) 153166
165
e ...... e Variable Obrechkoff Estimate e6 RungeKutto Estimate O 0 Exact Solution
1.00
0.75
0.50
0.25 o
0.00 o
0.25
0.50
0.75
V
1.00
2,00
0,00
4,00
6.00
8.00
10.00
Time
(The
exact solution is plotted using more points than ore shown)
1.00
0.75
I
"
"
!~
,
I,
o
,
t
~'
','
0.50
0.25 (3
0.00 (3_
ii'
,
,
,
i
!
I
I
i
i
i
i
! i
I t iIi?ititlil!l;ti
0.25
0.50
0,75
1.00
I
12.00
I
i
I
I
I
I
I
14.00 Time
[
I
I
I
I
i
i
I
16.00
Fig. 2. A comparison of variable Obrechkoff and classical fourstage RungeKutta estimates, for parameter a, using step size h = ~,' in the rigid rotating body problem [ p q r] T = [4 4t 2t] r, [a o b o c o do] r = [1 0 0 0] T.
166
D.P. Stapleton / Computer Physics Communications 98 (1996) 153166
=  ~2 otherwise, with the classical fourstage RungeKutta estimate and the (approximate by RungeKutta integration with a smaller step size) actual solution for parameter a. Both methods have fifth order local truncation error. The RungeKutta estimate, however, decays while the Obrechkoff estimate follows the solution's amplitude  the Obrechkoff plot is jagged at later times because the step size to period ratio is large, not because of an oscillating amplitude. Evidence for the superior stability of the Obrechkoff method was provided by RungeKutta runs (not shown) with smaller step sizes, in which estimates of the parameter a were initially much more accurate but eventually decayed.
5. Conclusions The implicit 1step variable Obrechkoff method (14) has been constructed to obtain the orientation in an inertial frame of a rotating rigid body for which the solution of Euler's equations is given. The theorem in Section 2 establishes the validity of the method in terms of stability, local truncation error, and maintaining a pure rotation. The method requires one evaluation of angular velocity components and one evaluation of angular acceleration components at each step. If the solution to Euler's equations is not given, (14), (24) and (25) can be used in an iterative procedure to attempt convergence simultaneously to the variable Obrechkoff estimate of the solution to Euler's equations and to the variable Obrechkoff estimate for the body orientation in an inertial frame. Eq. (11) indicates that implicit 1step variable Obrechkoff integrations with different variable coefficients than (12) are possible   with similar capabilities for absolute stability, local truncation error, and ability to maintain a pure rotation.
References [1] F.J. Regan, Reentry Vehicle Dynamics, AIAA Education Series (AIAA, New York, 1984) pp. 397408. [2] J.D. Lambert, Computational Methods in Ordinary Differential Equations (Wiley, reprinted by Arrowsmith, Bristol, 1979) p. 199.