# Solving Rotational Equations of Motion

I'm trying to integrate the rotational equations of motion

jTime = 1.0;
jMass = {{2.6666, 0, 0}, {0, 2.6666, 0}, {0, 0, 2.6666}};
jForce = (time^2*{-1.0, 4.0, -10.0}) + {2.0, -4.0, 10.0}
jAccel = jForce . Inverse[jMass]
jVel0 = {1, 1, 1};
jOri0 = {1, 0, 0, 0}
solOri =
NDSolve[{
jOri''[time] == jAccel,
jOri' == 0.5*Join[{0}, jVel0]*jOri0,
jOri == jOri0},jOri,
{time, 0, jTime}]


But I get the error

NDSolve::ndincd: "Initial conditions for derivatives of the function jOri[time] do not have consistent dimensions."

Well, this makes sense. The angular acceleration and velocity are 3D vectors while the orientation is a quaternion. Its also a little tricky to go between an angular velocity and an orientation, so we can't expect Mathematica to work that out by itself. My question is how to explain this to Mathematica?

EDIT: I tried doing everything as Euler angles

jForce = Cross[{1.0, 0.0, 0.0}, {0.0, -0.5*[Pi], 0.0}]
jRot = RotationMatrix[Pi - jOri[time][], {0, 0,    1}].RotationMatrix[jOri[time][], {1, 0, 0}].RotationMatrix[Pi - jOri[time]   [], {0, 0, 1}]
jAccel = jForce. Inverse[Transpose[jRot].jMass.jRot]
jOri0 = {0, 0, 0}
solOri = NDSolve[{jOri''[time] == jAccel, jOri' == jVel0, jOri == jOri0}, jOri, {time, 0, jTime}]


This was OK until I tried rotating the inertia tensor into the current frame. Now I have an error saying jOri[time][] doesn't exist. This seems contrary to all my initial conditions being 3 vectors

• If jOri is a quaternion, then jOri' will also be quaternion. You can't always mix notation that has different dimensions. Either change the orientation to cartesian directions, or change the equations of motion to quaternean form. – lalmei Feb 16 '14 at 12:07
• Normally I would convert to quaternion at the angular velocity stage (as amended above). However, I'm not sure how or if I can transform the accelation – RichardBruce Feb 16 '14 at 14:04
• while the quaternion for the angular velocity is right, I don't think you can multiply them that way. You need to use quaternion algebra. The dot product has a metric that I think is like a clifford algebra, I may be wrong. I will look into it. – lalmei Feb 17 '14 at 11:42

Let' s do this problem using just quaternions. First, let's see how we multiply quaternions together, remember that each component has those "imaginary" objects that multiply in a specific way, this can be replaced by just defining a new quarternion cross product that behaves the same way. We will keep the differential equations real so we will use the 4 dimensional representation :

$w\times q = \begin {pmatrix} w_ 0 & - w_ 1 & - w_ 2 & - w_ 3 \\ w_ 1 & w_ 0 & - w_ 3 & w_ 2 \\ w_ 2 & w_ 3 & w_ 0 & - w_ 1 \\ w_ 3 & - w_ 2 & w_ 1 & w_ 0 \end {pmatrix}\cdot q$

where $w = (w_ 0, w_ 1, w_ 2, w_ 3)$ and $q = (q_ 0, q_ 1, q_ 2, q_ 3)^T$, the vector components of the matrix ${M} _ {ij}$ is just $\epsilon_ {ijk} w_k$, just like a cross product, while the one that have components ${} _ 0$ are a combination of a asymmetric identity and the identity matrix, I wont go into too much detail as to how this comes about (maybe later, or look for 4 dimensional representation of Pauli matrices.)

In Mathematica we can simply write as

wquartpredot[w_] :=
Quiet@SparseArray[{{i_, i_} ->
w[], {i_, j_} /; i == 1 -> (-1)^i w[[j]], {i_, j_} /;
j == 1 -> (1)^j w[[i]],
{i_, j_}/;i != j -> (w[[2 ;;]].LeviCivitaTensor)[[j - 1, i - 1]]}, {4, 4}];


Now the $\times$ products are just wquartpredot[w].Transpose[q]

The equation of motion are also different.Since we are solving for the quaternions.

$$\dot q = \frac{1}{2} q \times \omega$$

as you wrote, but note that the product is defined as above.The second derivative now has two terms.

$$\ddot q = \frac{1}{2}\left ( \dot q\times\omega + q\times \dot \omega \right)$$

As long as the Inertia matrix is diagonal then

$$\dot\omega = \mathbf {I}^{-1} \tau$$

where $\tau$ is the torque. Now another quick note about dot products and $\times$ products.

Another thing we get from quaternions is conjugation.

$$\bar q = ( q_0,-q_1,-q_2,-q_3)$$

In Mathematica we implement this as

qconj[v_]:={{1,0,0,0},{0,-1,0,0},{0,0,-1,0},{0,0,0,-1}}.v;


Note that

$$q \times \bar q = ( q_0^2 +q_1^2+q_2^2+q_3^2, 0, 0, 0),$$

$$q \cdot \bar q = q_0^2 +q_1^2+q_2^2+q_3^2,$$

and more importantly for us

$$\bar q \times ( q \times v) = (\bar q \cdot q) \; v.$$

This allows to rewrite all the angular velocities in terms of the quaternions

$$\omega = 2 \bar q \times \dot q$$

We now need to solve the following equation, with its initial conditions,

$$\ddot q = \frac{1}{2}\left ( 2 \dot q\times ( \bar q \times \dot q) + q\times ( \mathbf {I}^{-1} \tau ) \right)$$

Assuming your jForce is the torque.

jTime = 1.0;
jMass = {{2.6666, 0, 0}, {0, 2.6666, 0}, {0, 0, 2.6666}};
jForce = (t^2*{-1.0, 4.0, -10.0}) + {2.0, -4.0, 10.0};
jAccel = Prepend[Inverse[jMass].jForce, 0];
JOri = {q0[t], q1[t], q2[t], q3[t]};
jVel0 = {0, 1, 1, 1};
jOri0 = {1, 0, 0, 0};


We write the equations above in mathematica

First the second derivative, quick and dirty:

doubledot=Table[D[D[JOri[[i]], t], t] == (1/2.*
(wquartpredot[JOri].Transpose[{jAccel}] +
wquartpredot[D[JOri, t]].(wquartpredot[qconj[JOri]].Transpose[{D[JOri,t]}]
)))[[i, 1]], {i, 1, 4}];


then the initial conditions

initialquart=MapThread[Equal, {{q0, q1, q2, q3}, jOri0}];

Derivative[q2], Derivative[q3]},
First@Transpose[wquartpredot[jOri0].Transpose[{jVel0}]]}];


Putting it all together, with a bit of massaging.

 sols=NDSolve[Flatten[Join[FullSimplify[doubledot],
{initialquart,initialdotquart}],1],{q0,q1,q2,q3},{t,0,jTime}];

{{q0->InterpolatingFunction[{{0.,1.}},<>],
q1->InterpolatingFunction[{{0.,1.}},<>],
q2->InterpolatingFunction[{{0.,1.}},<>],
q3->InterpolatingFunction[{{0.,1.}},<>]}}


Now at least gives a solution :)

• Looks good, fixed a couple of bugs. Roughly matches my own integrator, but runs forever if I try to rotate the inertia tensor to the current orientation, {q0,q1,q2,q3} – RichardBruce Feb 18 '14 at 14:02
• Thanks for catching the typos. If you are using a non-diagonal inertia tensor, I think you can't use the simple replacement for $\dot \omega$, there maybe extra terms proportional to $\omega$ I believe. I will look into later today. – lalmei Feb 18 '14 at 14:08
• Thanks, but I should be alright with your code now unless of course you just really enjoy the challenge :) – RichardBruce Feb 19 '14 at 1:26