# The Solution ofSO(3) through a Single Parameter ODE

Authors: Chein-Shan Liu
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

In many applications we need to solve an orthogonal transformation
tensorQ∈SO(3) from a tensorial equation Q˙=WQ under a given spin history W. In this paper, we address some interesting issues about this equation.
A general solution of Q is
obtained by transforming the governing equation into a new one in the space of ℝP3.
Then, we develop a novel method to solve Q in terms of
a single parameter, whose governing equation is a single nonlinear
ordinary differential equation (ODE).


---

## Body

## 1. Introduction

Among many classical Lie groups, the three-dimensional rotation groupSO(3) is the most widely used one. For its numerous engineering applications the development of simpler algorithms to calculate SO(3) under a large rotation has received a considerable attention in the literature. A comprehensive review on the spacecraft attitude was given by Shuster [1] and on the solid mechanics by Atluri and Cazzani [2]. A framework of minimal parameterizations of the rotation matrix was proposed by Bauchau and Trainelli [3].The purpose of searching for a suitable spin tensor, in a word, is to find a reference configuration with zero spin throughout the whole motion, such that the constitutive equation for a rate-type material under large deformation can be objectively integrated. To characterize this spin-free reference configuration/corotational frame, an orthogonal transformation tensorQ, connecting the spin-free and the fixed configurations due to the nonzero spin tensor denoted by W, satisfies the following tensorial differential equation:
(1)Q˙=WQ.
It does not lose any generality to assume that the initial condition of Q is an identity; that is, Q(0)=I3. Throughout this paper, a superimposed dot denotes the differential with respect to the current time t. Computational techniques were proposed by Rubinstein and Atluri [4] for integrating (1), which required a constant rate of rotation for each time step.It should be noted that the history ofQ can be represented by the histories of three Euler's angles θ, ϕ, and ψ as follows [5]:(2)Q=(cos⁡ψcos⁡ϕ-cos⁡θsinϕsinψcos⁡ψsinϕ+cos⁡θcos⁡ϕsinψsinψsinθ-sinψcos⁡ϕ-cos⁡θsinϕcos⁡ψ-sinψsinϕ+cos⁡θcos⁡ϕcos⁡ψcos⁡ψsinθsinθsinϕ-sinθcos⁡ϕcos⁡θ),and the corresponding differential equations are
(3)ω1=θ˙cos⁡ϕ+ψ˙sinθsinϕ,ω2=θ˙sinϕ-ψ˙sinθcos⁡ϕ,ω3=ϕ˙+ψ˙cos⁡θ.
Provided that the angular velocities ω1, ω2, and ω3 are given, the above nonlinear ordinary differential equations (ODEs) need to be integrated in a time-marching direction.For an effective representation of the rotation matrix, it has led to the development of numerous techniques in the last several decades, and the review of the properties, advantages, and shortcomings of these parameterization techniques can be found in Ibrahimbegovic [6], Borri et al. [7], and Bauchau and Trainelli [3]. To represent the three-dimensional rotation, usually the number of parameters is three, like the Euler parameters, the Rodrigues parameters, and the modified Rodrigues parameters. However, these representations contain certain singularities, and their governing equations are highly nonlinear in nature. The procedures for finding the solutions of rotation matrix involving these nonlinear ODEs systems are usually very complicated.It is known that the spatial orientationQ∈SO(3) of a rigid body rotation can be expressed in terms of the unit quaternion [8]:(4)Q=(q02+q12-q22-q322(q1q2+q0q3)2(q1q3-q0q2)2(q1q2-q0q3)q02-q12+q22-q322(q2q3+q0q1)2(q1q3+q0q2)2(q2q3-q0q1)q02-q12-q22+q32).These parameters are obtained by using the stereographic projection of
(5)𝕊3∶={q=(q,q0)∈ℝ4∣∥q∥2+q02=1}
onto ℝ3 by a two-fold covering; see, for example, Goldstein [5]. In the above, ∥q∥ denotes the Euclidean norm of q∈ℝ3. Liu [8] has presented a four-dimensional Lie-algebra representation of the quaternion formulation of SO(3).It is known thatSU(2) is diffeomorphic to the three-dimensional sphere 𝕊2 and SO(3) is diffeomorphic to the quotient space of the three-dimensional sphere by the antipodal equivalence, hence diffeomorphic to the three-dimensional projective space [9]. According to [10] we can define the real projective space as follows. Let ℝPn be the set of all straight lines through the origin in ℝn. a and b∈ℝn represent the same line if and only if there is a nonzero constant λ∈ℝ such that a=λb, which constitutes an equivalence class denoted by [a]={b∣b~a} with b~a meaning that a=λb for some λ≠0. Note that ℝPn=(ℝn-{0})/~. The coordinates of any a∈ℝn such that b=[a] are called homogeneous coordinates for b.In this paper, a simpler solution method ofQ is proposed by expressing the orientation equations in terms of local coordinates, yielding a scalar differential equation which can avoid the singularity. For a specified level of accuracy in numerical integration, a scalar equation requires less CPU time than an equivalent transcendental set of ODEs as shown in (3). To interpret the results of integration, the time evolution of orientation is presented as a curve with a single parameter in the three-dimensional topological space ℝP3. The local coordinate is an example of a globally defined nonsingular parameterization of rotations, which is suitable for treating the computations of large rotations.


## 2. A Decomposition ofQ

We denote the spin matrix by(6)W=(0-ω3ω2ω30-ω1-ω2ω10),
and the corresponding angular velocity vector is
(7)ω=(ω1,ω2,ω3)T,
whose magnitude is denoted by
(8)∥ω∥∶=(ω12+ω22+ω32)1/2.
Also the instantaneous spin axis in the three-dimensional space is denoted by
(9)(η1,η2,η3)T:=1∥ω∥(ω1,ω2,ω3)T.
When the spin axis is fixed, it can be viewed as a two-dimensional (2D) spin since the rotation only occurs on the plane which is perpendicular to a fixed axis. While the spin axis is varying with time, it is a three-dimensional (3D) spin.In what follows, we present a novel method to explore the general solution ofQ. Since Q is orthogonal, it belongs to the special orthogonal group with dimensions three; that is, Q∈SO(3). Although (1) can be defined by nine simultaneous ordinary differential equations (ODEs), only three of them are independent. In geometry, Q, an element of SO(3), represents a certain 3D algebraic surface in a real space of nine dimensions. It is unwise to find the analytical solution of (1) by solving these simultaneous ODEs. Here, the problem is solved by a judicious consideration based on a novel technique. For the sake of convenience, let us define a matrix operator F which applies to W and has the following form:
(10)F(W):=I3+sinw∥ω∥W+(1-cos⁡w)∥ω∥2W2,
where w(t):=∫0t∥ω(τ)∥dτ.We consider a subset of (6), by defining the following 2D spin matrix:
(11)W1=(00000-ω10ω10),
from which we have
(12)F(W1)=(1000cos⁡ω¯1-sinω¯10sinω¯1cos⁡ω¯1),
where
(13)ω¯1(t):=∫0t‍ω1(τ)dτ.It is cunning to presume that the solution ofQ can be decomposed into
(14)Q=F(W1)Q1,
with Q1 an unknown matrix belonging to SO(3). Substituting it into (1) leads to
(15)Q˙1=AQ1,Q1(0)=I3,
where
(16)A:=F(W1)T[WF(W1)-F˙(W1)]=(0-u˙1-u˙2u˙100u˙200)
is a skew-symmetric matrix with
(17)u˙1:=ω3cos⁡ω¯1-ω2sinω¯1,u˙2:=-ω2cos⁡ω¯1-ω3sinω¯1.The decomposition in (14) leads to a simpler spin matrix A for Q1 in (16), with only two independent inputs u˙1 and u˙2. For a given angular velocity (ω1(t),ω2(t),ω3(t)), it is easy to find the matrix F by (13) and (12). However, in order to obtain Q we still require to find Q1. In this paper, an analytic procedure will be developed to solve this problem for arbitrary inputs u˙1 and u˙2 generated from the angular velocity (ω1(t),ω2(t),ω3(t)).


## 3. A Projective Transformation

The system of ODEs deduced fromA in (16) can be written as
(18)X˙0=-u˙1X1-u˙2X2,X˙1=u˙1X0,X˙2=u˙2X0.
The initial values of X0, X1, and X2 are assumed to be X0(0), X1(0), and X2(0), respectively. So the determination of Q1(t) is now equivalent to searching a general solution of (18); that is,
(19)X(t)=Q1(t)X0,
where Q1(0)=I3, X0=X(0), and
(20)X(t):=(X0(t)X1(t)X2(t)).Let(21)x1:=X1X0,x2:=X2X0
be the homogeneous coordinates of ℝP3. Then, the use of (18) implies
(22)X˙0X0=-x·u˙,(23)ddt(X0x)=X0u˙,
where
(24)x:=(x1x2),u˙:=(u˙1u˙2)
are the output and input of (22) and (23), respectively.The inner product of (23) with x and the use of (22) render
(25)x·x˙=-X˙0X0(∥x∥2+1).
Integrating (25) leads to
(26)∥x(t)∥2=∥X0∥2X02(t)-1.
By (21), it is equivalent to ∥X(t)∥=∥X0∥; that is, the length of the vector X is preserved under the action of SO(3) group. Obviously, X0 cannot be a zero vector; otherwise, X(t) will be a zero vector for all t>0.By eliminatingX0, (22) and (23) can be combined into a nonlinear differential equations system for x:
(27)x˙-(u˙·x)x=u˙.
The transformation made in this section projects the three-dimensional vector (X0,X1,X2)T∈𝕊∥X0∥2, where 𝕊∥X0∥2 means a three-dimensional sphere with a constant radius ∥X0∥, into a two-dimensional vector (x1,x2)T in the topological space ℝP3, which is correlated intimately with the two independent inputs of u˙.


## 4. The Main Results

### 4.1. Two Theorems

In this section we are going to prove two main theorems.Theorem 1.
The solution ofx governed by (27) with an initial condition x(0)=x0 can be explicitly expressed in terms of a single variable z¯:
(28)x=z(z¯)cos⁡z¯-x0·v˙sinz¯,(29)z(z¯)=x0+[sinz¯+x0·v˙(cos⁡z¯-1)]v˙,
where v is a constant unit vector (with ∥v˙∥=1), and z¯ is governed by a nonlinear ODE:
(30)z¯˙=(v˙1u˙2-v˙2u˙1)[∥z∥2+(cos⁡z¯-x0·v˙sinz¯)2](sinz¯+x0·v˙cos⁡z¯)[x2(0)v˙1-x1(0)v˙2]+z2(z¯)u˙1-z1(z¯)u˙2x2(0)v˙1-x1(0)v˙2,
under the initial condition z¯(0)=0.Proof.
The proof of this theorem is quite lengthy, and we divide it into five parts.
(A) Mixing the Input and Output. Consider the following transformations of variables:
(31)x˙=w˙+v˙,(32)u˙=cw˙+dv˙,
where v˙ is aconstant vector with norm ∥v˙∥=1 to be given, and the vector w˙ and the other two scalars c and d are allowed to betime-varying. We will determine v˙, w˙, c, and d below, under the assumptions c≠d and c≠1.
Substituting (31) and its integral into (23) we can obtain
(33)w˙+X˙0(1-c)X0w=d-11-cv˙-X˙0(1-c)X0v,
where v(t):=tv˙.
Upon defining(34)y˙y=X˙0(1-c)X0,
equation (33) becomes
(35)w˙+y˙yw=d-11-cv˙-y˙yv,
the solution of which is
(36)w(t)=y(0)y(t)w0+∫0t[d(τ)-1]y(τ)[1-c(τ)]y(t)v˙dτ-∫0ty˙(τ)y(t)v(τ)dτ,
where w0=w(0)=x0. The last term can be integrated by parts, leading to
(37)w(t)=-v(t)+y(0)y(t)x0+v˙∫0t‍z(τ)y(τ)y(t)dτ,
where
(38)z:=d-c1-c
is a time function. Under the conditions c≠1 and c≠d, z is a well-defined nonzero function. It is remarkable that (37) expresses w in terms of a constant unit vector v˙.
(B) Governing Equations of
w
and
y. From (25) and (34) it follows that
(39)x·x˙∥x∥2+1=(c-1)y˙y.
The inner product of w+v and (35) is
(40)w˙·(w+v)+y˙y(w+v)·(w+v)=d-11-c(w+v)·v˙,
and by (31) we have
(41)w˙·(w+v)+v˙·(w+v)=x·x˙,(w+v)·(w+v)=∥x∥2.
Thus (40) can be changed to
(42)x·x˙+y˙y∥x∥2=z(w+v)·v˙,
which upon using (39) becomes
(43)(c∥x∥2+c-1)y˙y=z(w+v)·v˙.
Without losing any generality we may select c as
(44)c=1-b∥x∥2+1,
where b>0 is a time function to be determined; hence, (43) becomes
(45)y˙y=-zb(w+v)·v˙.
Equations (35) and (45) are composed as the governing equations system for (w,y), with v being the input.
(C) Explicit Form of
y. Noting (37), (45) changes to
(46)y˙(t)=-y(0)z(t)b(t)x0·v˙-1b(t)∫0t‍z(t)z(τ)y(τ)dτ.
Define
(47)t′:=∫0t‍1b(τ)dτ,
and the relation between t′ and t is one-to-one, since b>0. Now, y is viewed as a function of t′, such that
(48)dy(t′)dt′=-y(0)z(t′)x0·v˙-∫0t′‍z(t′)z(s)y(s)ds,
by (46) and (47).
From (48) we have
(49)y′(z¯(t′))=-y(0)x0·v˙-∫0z¯(t′)‍y(z¯(s))dz¯(s),
where y′ denotes the differential with respect to z¯, which is defined by
(50)z¯(t′)≔∫0t′‍z(s)ds.
Taking the differential of (49) with respect to z¯ again, we can obtain
(51)y′′=-y.
The solution of y is
(52)y=cos⁡z¯-x0·v˙sinz¯,
where y(0)=1 and y′(0)=-x0·v˙ are imposed. It is interesting that we have a closed-form solution of y in terms of z¯.
(D) Explicit Form of
x
. Now, substituting (52) into (37) and integrating the resultant, the explicit form of x can be obtained as follows:
(53)x=w+v=zy,
where
(54)z:=x0+sinz¯v˙+x0·v˙(cos⁡z¯-1)v˙.
Ifz¯ can be solved, the solution of x is obtained. Defining the vector z as given in (29) and substituting (52) for y into the above equation, we obtain (28). The square norm of x is given by
(55)∥x∥2=∥z∥2y2.
(E) The Governing Equation of
z
¯. It can be seen that the single parameter of variable z¯ plays the major role to express the solutions of y and x above. The issue to find z¯ is very important as being given below. Using (35), (38), and (53) we have
(56)w˙=-y˙y2z+(z-1)v˙.
Substituting (56) for w˙ into (32) and using (38), one has
(57)u˙=-cy˙y2z+zv˙.
From (44) and (55) the term c reads as
(58)c=(1-b)y2∥z∥2+y2,
which together with
(59)z=bz¯˙,
a result deduced from (50) and (47), and y˙=y′z¯˙ being substituted into (57), renders
(60)u˙=(b-1)y′z¯˙z∥z∥2+y2+bz¯˙v˙.
In component form we have
(61)u˙1=(b-1)y′z¯˙z1∥z∥2+y2+bz¯˙v˙1,(62)u˙2=(b-1)y′z¯˙z2∥z∥2+y2+bz¯˙v˙2.
The above two equations can be used to solveb and z¯˙. Eliminating z¯˙ from the above two equations we can obtain
(63)b(t)=F1(z¯)u˙2-F2(z¯)u˙1v˙1u˙2-v˙2u˙1+F1(z¯)u˙2-F2(z¯)u˙1,
where
(64)F1(z¯):=y′z1∥z∥2+y2,F2(z¯):=y′z2∥z∥2+y2.
It can be seen that b is a function of z¯ and t, and the latter is induced by the inputs u˙1 and u˙2. Multiplying (62) by F1 and (61) by F1 and then subtracting them we obtain
(65)z¯˙=F1(z¯)u˙2-F2(z¯)u˙1b[F1(z¯)v˙2-F2(z¯)v˙1].
After substituting (63) for b into the above equation we can derive
(66)z¯˙=[v˙1+F1(z¯)]u˙2-[v˙2+F2(z¯)]u˙1F1(z¯)v˙2-F2(z¯)v˙1.
It is a first order ODE for z¯ under the initial condition z¯(0)=0, the integration of which gives z¯(t). With the aid of (64), (52), and (29) and through some manipulations the above equation leads to (30). This ends the proof of Theorem 1.Theorem 2.
The solutions ofX are represented by
(67)X0(t)=(cos⁡z¯-x0·v˙sinz¯)∥X0∥∥z∥2+(cos⁡z¯-x0·v˙sinz¯)2,(68)X1(t)=z1(z¯)∥X0∥∥z∥2+(cos⁡z¯-x0·v˙sinz¯)2,(69)X2(t)=z2(z¯)∥X0∥∥z∥2+(cos⁡z¯-x0·v˙sinz¯)2.Proof.
From (63), (52), and (65) the history of b can be obtained. Thus, from (58) and (59) the histories of c and z can be calculated, respectively, whereas through (26) and (55) the history of X0 can be evaluated as follows:
(70)X0(t)=y∥X0∥∥z∥2+y2.
Inserting (52) for y into the above equation we can derive (67). The last two components X1 and X2 of X can be obtained explicitly via (21), (67), and (28), from which we can derive (68) and (69). This ends the proof of Theorem 2.From (67)–(69) it is obvious that ∥X(t)∥=∥X0∥; that is, (X0(t),X1(t),X2(t))T∈𝕊∥X0∥2, where 𝕊∥X0∥2 means a three-dimensional sphere with a constant radius ∥X0∥. Therefore, the above mappings belong to SO(3), which preserves the invariant length of X. The above result is significant upon recalling the number of parameters used in the Euler's angels is three and the governing equations are three nonlinear ODEs. Here, we only need a single parameter z¯ and a nonlinear ODE in (66).


### 4.2. The Choice ofv

If we choose(71)v˙1=u˙1(0)u˙12(0)+u˙22(0),v˙2=u˙2(0)u˙12(0)+u˙22(0),
we have b(0)=1 by (63). In (29) we expand the two-dimensional vector z in the plane by two oblique coordinates x0 and v˙, which are subjected to the constraint that they are not parallel; that is, x2(0)v˙1-x1(0)v˙2≠0; otherwise, the denominators in (30) are zero, which would lead to an undefined differential equation for z¯. In the whole process b(t)>0 should hold. Upon b(t)<e0, where 0<e0<1 is a given constant, at some time instant, say t0, the numerical integration process is restarted with new values of v˙1 and v˙2 given by
(72)v˙1=u˙1(t0)u˙12(t0)+u˙22(t0),v˙2=u˙2(t0)u˙12(t0)+u˙22(t0),
and at the same time we set z¯(t0)=0, x1(t0)=X1(t0)/X0(t0), and x2(t0)=X2(t0)/X0(t0).


### 4.3. The Computation ofQ

First we calculateQ1 by the above method. Select three independent initial values of X0; for example,
(73)(X01(0)X02(0)X03(0)X11(0)X12(0)X13(0)X21(0)X22(0)X23(0))=(111110101).
The corresponding solutions are denoted by
(74)(X01(t)X02(t)X03(t)X11(t)X12(t)X13(t)X21(t)X22(t)X23(t)),
and from (19) we obtain(75)Q1(t)=(X01(t)X02(t)X03(t)X11(t)X12(t)X13(t)X21(t)X22(t)X23(t))(111110101)-1=(X02(t)+X03(t)-X01(t)X01(t)-X03(t)X01(t)-X02(t)X12(t)+X13(t)-X11(t)X11(t)-X13(t)X11(t)-X12(t)X22(t)+X23(t)-X21(t)X21(t)-X23(t)X21(t)-X22(t))Then, inserting the above equation for Q1(t) and (12) for F into (14) we can obtain Q(t):(76)Q(t)=(1000cos⁡ω¯1-sinω¯10sinω¯1cos⁡ω¯1)(X02(t)+X03(t)-X01(t)X01(t)-X03(t)X01(t)-X02(t)X12(t)+X13(t)-X11(t)X11(t)-X13(t)X11(t)-X12(t)X22(t)+X23(t)-X21(t)X21(t)-X23(t)X21(t)-X22(t)).


### 4.4. Numerical Tests

In order to give a criterion to assess our numerical method we first derive a closed-form solution ofQ1 in the appendix under the angular velocities ω1=Ω-ω, ω2=-sinΩt, and ω3=cos⁡Ωt, where Ω and ω are parameters of angular frequencies.We calculate(X0(t),X1(t),X2(t)) by the above single-parameter method with e0=0.999. The initial value is (X0(0),X1(0),X2(0))=(5,2,10) and we fix ω=2, Ω=3. We apply the Euler method to integrate (30) by using a stepsize h=0.0001. The errors between exact solutions and numerical solutions are plotted in Figure 1, whose maximum errors are 1.55×10-4, 5.91×10-4, and 2.1×10-4, respectively.Error between the exact and the numerical solution provided by the single-parameter method.
(a)
(b)
(c)Now, let us turn to the case of a large rotation in Figure2 up to t=10, where ω=10 and Ω=5 were used. We apply the Euler method to integrate (30) by using a stepsize h=0.001 and e0=0.96. Then the method in Section 4.3 is used to compute Q, whose componential errors are shown in Figures 2(a)–2(i), where the maximum error is smaller than 10-3. The error of orthogonality is defined as ∥QTQ-I3∥ with Q calculated by the numerical method. From Figure 2(j) it can be seen that the present numerical method can preserve orthogonality almost exactly.Component-wise errors of rotation matrix (a)–(i) and error of orthogonality (j) produced by the single-parameter method.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)


## 4.1. Two Theorems

In this section we are going to prove two main theorems.Theorem 1.
The solution ofx governed by (27) with an initial condition x(0)=x0 can be explicitly expressed in terms of a single variable z¯:
(28)x=z(z¯)cos⁡z¯-x0·v˙sinz¯,(29)z(z¯)=x0+[sinz¯+x0·v˙(cos⁡z¯-1)]v˙,
where v is a constant unit vector (with ∥v˙∥=1), and z¯ is governed by a nonlinear ODE:
(30)z¯˙=(v˙1u˙2-v˙2u˙1)[∥z∥2+(cos⁡z¯-x0·v˙sinz¯)2](sinz¯+x0·v˙cos⁡z¯)[x2(0)v˙1-x1(0)v˙2]+z2(z¯)u˙1-z1(z¯)u˙2x2(0)v˙1-x1(0)v˙2,
under the initial condition z¯(0)=0.Proof.
The proof of this theorem is quite lengthy, and we divide it into five parts.
(A) Mixing the Input and Output. Consider the following transformations of variables:
(31)x˙=w˙+v˙,(32)u˙=cw˙+dv˙,
where v˙ is aconstant vector with norm ∥v˙∥=1 to be given, and the vector w˙ and the other two scalars c and d are allowed to betime-varying. We will determine v˙, w˙, c, and d below, under the assumptions c≠d and c≠1.
Substituting (31) and its integral into (23) we can obtain
(33)w˙+X˙0(1-c)X0w=d-11-cv˙-X˙0(1-c)X0v,
where v(t):=tv˙.
Upon defining(34)y˙y=X˙0(1-c)X0,
equation (33) becomes
(35)w˙+y˙yw=d-11-cv˙-y˙yv,
the solution of which is
(36)w(t)=y(0)y(t)w0+∫0t[d(τ)-1]y(τ)[1-c(τ)]y(t)v˙dτ-∫0ty˙(τ)y(t)v(τ)dτ,
where w0=w(0)=x0. The last term can be integrated by parts, leading to
(37)w(t)=-v(t)+y(0)y(t)x0+v˙∫0t‍z(τ)y(τ)y(t)dτ,
where
(38)z:=d-c1-c
is a time function. Under the conditions c≠1 and c≠d, z is a well-defined nonzero function. It is remarkable that (37) expresses w in terms of a constant unit vector v˙.
(B) Governing Equations of
w
and
y. From (25) and (34) it follows that
(39)x·x˙∥x∥2+1=(c-1)y˙y.
The inner product of w+v and (35) is
(40)w˙·(w+v)+y˙y(w+v)·(w+v)=d-11-c(w+v)·v˙,
and by (31) we have
(41)w˙·(w+v)+v˙·(w+v)=x·x˙,(w+v)·(w+v)=∥x∥2.
Thus (40) can be changed to
(42)x·x˙+y˙y∥x∥2=z(w+v)·v˙,
which upon using (39) becomes
(43)(c∥x∥2+c-1)y˙y=z(w+v)·v˙.
Without losing any generality we may select c as
(44)c=1-b∥x∥2+1,
where b>0 is a time function to be determined; hence, (43) becomes
(45)y˙y=-zb(w+v)·v˙.
Equations (35) and (45) are composed as the governing equations system for (w,y), with v being the input.
(C) Explicit Form of
y. Noting (37), (45) changes to
(46)y˙(t)=-y(0)z(t)b(t)x0·v˙-1b(t)∫0t‍z(t)z(τ)y(τ)dτ.
Define
(47)t′:=∫0t‍1b(τ)dτ,
and the relation between t′ and t is one-to-one, since b>0. Now, y is viewed as a function of t′, such that
(48)dy(t′)dt′=-y(0)z(t′)x0·v˙-∫0t′‍z(t′)z(s)y(s)ds,
by (46) and (47).
From (48) we have
(49)y′(z¯(t′))=-y(0)x0·v˙-∫0z¯(t′)‍y(z¯(s))dz¯(s),
where y′ denotes the differential with respect to z¯, which is defined by
(50)z¯(t′)≔∫0t′‍z(s)ds.
Taking the differential of (49) with respect to z¯ again, we can obtain
(51)y′′=-y.
The solution of y is
(52)y=cos⁡z¯-x0·v˙sinz¯,
where y(0)=1 and y′(0)=-x0·v˙ are imposed. It is interesting that we have a closed-form solution of y in terms of z¯.
(D) Explicit Form of
x
. Now, substituting (52) into (37) and integrating the resultant, the explicit form of x can be obtained as follows:
(53)x=w+v=zy,
where
(54)z:=x0+sinz¯v˙+x0·v˙(cos⁡z¯-1)v˙.
Ifz¯ can be solved, the solution of x is obtained. Defining the vector z as given in (29) and substituting (52) for y into the above equation, we obtain (28). The square norm of x is given by
(55)∥x∥2=∥z∥2y2.
(E) The Governing Equation of
z
¯. It can be seen that the single parameter of variable z¯ plays the major role to express the solutions of y and x above. The issue to find z¯ is very important as being given below. Using (35), (38), and (53) we have
(56)w˙=-y˙y2z+(z-1)v˙.
Substituting (56) for w˙ into (32) and using (38), one has
(57)u˙=-cy˙y2z+zv˙.
From (44) and (55) the term c reads as
(58)c=(1-b)y2∥z∥2+y2,
which together with
(59)z=bz¯˙,
a result deduced from (50) and (47), and y˙=y′z¯˙ being substituted into (57), renders
(60)u˙=(b-1)y′z¯˙z∥z∥2+y2+bz¯˙v˙.
In component form we have
(61)u˙1=(b-1)y′z¯˙z1∥z∥2+y2+bz¯˙v˙1,(62)u˙2=(b-1)y′z¯˙z2∥z∥2+y2+bz¯˙v˙2.
The above two equations can be used to solveb and z¯˙. Eliminating z¯˙ from the above two equations we can obtain
(63)b(t)=F1(z¯)u˙2-F2(z¯)u˙1v˙1u˙2-v˙2u˙1+F1(z¯)u˙2-F2(z¯)u˙1,
where
(64)F1(z¯):=y′z1∥z∥2+y2,F2(z¯):=y′z2∥z∥2+y2.
It can be seen that b is a function of z¯ and t, and the latter is induced by the inputs u˙1 and u˙2. Multiplying (62) by F1 and (61) by F1 and then subtracting them we obtain
(65)z¯˙=F1(z¯)u˙2-F2(z¯)u˙1b[F1(z¯)v˙2-F2(z¯)v˙1].
After substituting (63) for b into the above equation we can derive
(66)z¯˙=[v˙1+F1(z¯)]u˙2-[v˙2+F2(z¯)]u˙1F1(z¯)v˙2-F2(z¯)v˙1.
It is a first order ODE for z¯ under the initial condition z¯(0)=0, the integration of which gives z¯(t). With the aid of (64), (52), and (29) and through some manipulations the above equation leads to (30). This ends the proof of Theorem 1.Theorem 2.
The solutions ofX are represented by
(67)X0(t)=(cos⁡z¯-x0·v˙sinz¯)∥X0∥∥z∥2+(cos⁡z¯-x0·v˙sinz¯)2,(68)X1(t)=z1(z¯)∥X0∥∥z∥2+(cos⁡z¯-x0·v˙sinz¯)2,(69)X2(t)=z2(z¯)∥X0∥∥z∥2+(cos⁡z¯-x0·v˙sinz¯)2.Proof.
From (63), (52), and (65) the history of b can be obtained. Thus, from (58) and (59) the histories of c and z can be calculated, respectively, whereas through (26) and (55) the history of X0 can be evaluated as follows:
(70)X0(t)=y∥X0∥∥z∥2+y2.
Inserting (52) for y into the above equation we can derive (67). The last two components X1 and X2 of X can be obtained explicitly via (21), (67), and (28), from which we can derive (68) and (69). This ends the proof of Theorem 2.From (67)–(69) it is obvious that ∥X(t)∥=∥X0∥; that is, (X0(t),X1(t),X2(t))T∈𝕊∥X0∥2, where 𝕊∥X0∥2 means a three-dimensional sphere with a constant radius ∥X0∥. Therefore, the above mappings belong to SO(3), which preserves the invariant length of X. The above result is significant upon recalling the number of parameters used in the Euler's angels is three and the governing equations are three nonlinear ODEs. Here, we only need a single parameter z¯ and a nonlinear ODE in (66).


## 4.2. The Choice ofv

If we choose(71)v˙1=u˙1(0)u˙12(0)+u˙22(0),v˙2=u˙2(0)u˙12(0)+u˙22(0),
we have b(0)=1 by (63). In (29) we expand the two-dimensional vector z in the plane by two oblique coordinates x0 and v˙, which are subjected to the constraint that they are not parallel; that is, x2(0)v˙1-x1(0)v˙2≠0; otherwise, the denominators in (30) are zero, which would lead to an undefined differential equation for z¯. In the whole process b(t)>0 should hold. Upon b(t)<e0, where 0<e0<1 is a given constant, at some time instant, say t0, the numerical integration process is restarted with new values of v˙1 and v˙2 given by
(72)v˙1=u˙1(t0)u˙12(t0)+u˙22(t0),v˙2=u˙2(t0)u˙12(t0)+u˙22(t0),
and at the same time we set z¯(t0)=0, x1(t0)=X1(t0)/X0(t0), and x2(t0)=X2(t0)/X0(t0).


## 4.3. The Computation ofQ

First we calculateQ1 by the above method. Select three independent initial values of X0; for example,
(73)(X01(0)X02(0)X03(0)X11(0)X12(0)X13(0)X21(0)X22(0)X23(0))=(111110101).
The corresponding solutions are denoted by
(74)(X01(t)X02(t)X03(t)X11(t)X12(t)X13(t)X21(t)X22(t)X23(t)),
and from (19) we obtain(75)Q1(t)=(X01(t)X02(t)X03(t)X11(t)X12(t)X13(t)X21(t)X22(t)X23(t))(111110101)-1=(X02(t)+X03(t)-X01(t)X01(t)-X03(t)X01(t)-X02(t)X12(t)+X13(t)-X11(t)X11(t)-X13(t)X11(t)-X12(t)X22(t)+X23(t)-X21(t)X21(t)-X23(t)X21(t)-X22(t))Then, inserting the above equation for Q1(t) and (12) for F into (14) we can obtain Q(t):(76)Q(t)=(1000cos⁡ω¯1-sinω¯10sinω¯1cos⁡ω¯1)(X02(t)+X03(t)-X01(t)X01(t)-X03(t)X01(t)-X02(t)X12(t)+X13(t)-X11(t)X11(t)-X13(t)X11(t)-X12(t)X22(t)+X23(t)-X21(t)X21(t)-X23(t)X21(t)-X22(t)).


## 4.4. Numerical Tests

In order to give a criterion to assess our numerical method we first derive a closed-form solution ofQ1 in the appendix under the angular velocities ω1=Ω-ω, ω2=-sinΩt, and ω3=cos⁡Ωt, where Ω and ω are parameters of angular frequencies.We calculate(X0(t),X1(t),X2(t)) by the above single-parameter method with e0=0.999. The initial value is (X0(0),X1(0),X2(0))=(5,2,10) and we fix ω=2, Ω=3. We apply the Euler method to integrate (30) by using a stepsize h=0.0001. The errors between exact solutions and numerical solutions are plotted in Figure 1, whose maximum errors are 1.55×10-4, 5.91×10-4, and 2.1×10-4, respectively.Error between the exact and the numerical solution provided by the single-parameter method.
(a)
(b)
(c)Now, let us turn to the case of a large rotation in Figure2 up to t=10, where ω=10 and Ω=5 were used. We apply the Euler method to integrate (30) by using a stepsize h=0.001 and e0=0.96. Then the method in Section 4.3 is used to compute Q, whose componential errors are shown in Figures 2(a)–2(i), where the maximum error is smaller than 10-3. The error of orthogonality is defined as ∥QTQ-I3∥ with Q calculated by the numerical method. From Figure 2(j) it can be seen that the present numerical method can preserve orthogonality almost exactly.Component-wise errors of rotation matrix (a)–(i) and error of orthogonality (j) produced by the single-parameter method.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)


## 5. Conclusions

Upon comparing with some different representations of the rotation groupSO(3), including the Euler's angles representation, the Rodrigues parameters representation, and the modified Rodrigues parameters representation, we succeeded to develop a simpler mathematical procedure to find an analytical solution of Q through a single parameter, where we just need to solve a single nonlinear ODE. To interpret the results of the integration, the time evolution of orientation is presented as a curve with a single parameter in the topological space ℝP3. The new local coordinate is a globally defined nonsingular parameterization of rotations suitable for general solutions of large rotations.



---

*Source: 671921-2014-04-30.xml*