How can I apply Euler's Method to predict a point in time rotating around multiple axis'

I am xposting this from my original stackoverflow question where I was presented with a coding challenge that I have been able to narrow down extensively and I think it lies with Euler's Method.

Here's my initial post with my previous assumptions.

It has been a while since I took differential equations and linear algebra so I am struggling to remember/keep up with the annotations of this and I am hoping someone can help me out.

\begin{aligned}I_{1}{\dot {\omega }}_{{1}}+(I_{3}-I_{2})\omega _{2}\omega _{3}&=M_{{1}}\\I_{2}{\dot {\omega }}_{{2}}+(I_{1}-I_{3})\omega _{3}\omega _{1}&=M_{{2}}\\I_{3}{\dot {\omega }}_{{3}}+(I_{2}-I_{1})\omega _{1}\omega _{2}&=M_{{3}}\end{aligned} The information I have:

Angular Velocity of the asteroid from the center of mass given in radians/second [wX,wY,wZ].

Current (Instantaneous) Position of the point on the asteroid at time(t)=0. This is with respect to the center of mass of the asteroid. Just a vector pointing away from the center of mass to the point I am wanting to hit.

A rough approximation of "t", such as 5.8 seconds. I am deriving the time (t) from the distance of the shooter to the center of mass of the asteroid. I may be able to get a better approximation of t if I can figure out the above.

The Problem:

I am developing a space shooter that works in 3D space. I am wanting the AI to be able to fire a projectile at a point on a tumbling asteroid so that when the projectile hits the asteroid, it's on that point.

If the asteroid is not rotating, it's a pretty simple kinematics equation with no differentials needed. However in some cases the asteroid is rotating around the X, Y, and Z axis. I am wanting to find a way to describe the motion of the asteroid and project where that point will be at a given time(t).

I have tried using 4x4 rotation/translation matricies for the above. This works if the asteroid is rotating around a single axis (say the Y axis), but the second it is rotating around more than 1 axis it looks like the rotational matricies are off (almost looks like its the opposite direction of what it should be)

The Information I am wanting to get

I am wanting to have a simple model where I can essentially plug in these values (angular velocity, current rotation, current position of the point, and time(t)) to get a point in space where that point will be at time(t).

I am struggling to teach myself this method. I'm sure it is more simple than I am thinking it is, I guess I'm struggling to understand what numbers go where and how.

If the initial values were (arbitrarily):

Angular Velocity = (1,-2,3) (in radians/second)

Current Point Position (with respect to center of mass of asteroid) = (5,4,6)

time(t)=7 seconds

How could I use this method to predict the location of the point at time (t) with respect to the asteroid's center of mass? Is it as simple as plugging in these values into the equations? Maybe that's where I'm overthinking this. As of now I am not accounting for angular acceleration, would it be exponentially more difficult if I did? (Say the asteroid had a weird thrusting booster on it)

Thank you so much for any insight you might help me with,

Chris

• You can decouple the linear and rotational motion. I assume there are no forces acting on the asteroid until it's hit? In that case the linear motion will be a straight line. I believe there is an analytic solution to the torque free Euler's equation of rotation, but I am not sure if it will also give an analytic solution for the rotation matrix, so you might have to do numerical ode solving, unless you are rotating around only one of the principal axis. – fibonatic Jun 2 '18 at 11:02
• @Chris I hope you have seen the answer to your question posted below. – Futurologist Jun 6 '18 at 2:13
• @Futurologist I will have time tomorrow to hopefully apply this. This is perfect, thank you so much. I will absolutely post the code I derive from what you have given me. Thanks again – Chris Jun 7 '18 at 2:23

This came out a bit long, so first you may want to scroll down to the algorithm, skipping the explanations and the justifications. I have written a discrete dynamical model of the case when the asteroid rotates freely in space, twisting and turning, with precessing angular velocity (so there is angular acceleration), without any thrusting boosters. As you will see, the dynamics is non-trivial and interesting enough.

I can suggest a sketch but you may have to try to figure out the details. I will use standard $$3$$ dimensional vectors, orthogonal and skew symmetric matrices but one can do the absolutely the same thing using quaternions. I am not sure which one is computationally cheaper or more convenient but maybe there is not that much difference. Ether way, the philosophy is always the same and conceptually everything looks the same regardless of the usage of quaternions or $$3$$D vectors and matrices.

We will be working a lot with the set of all orthogonal matrices (also known as the orthogonal group) $$\text{SO}(3) = \Big\{ \, U \in \mathbb{R}^{3 \times 3} \,\, \big{|} \,\, U\, U^T = U^T U = I \, \Big\}$$ where $$U^T$$ is the transpose matrix of $$U$$. Furthermore, we will be working with the set of all skew-symmetric matrices (also known as the Lie algebra of $$\text{SO}(3)$$) $$\text{so}(3) = \Big\{ \, A \in \mathbb{R}^{3 \times 3} \,\, \big{|} \,\, A^T = -\, A \, \Big\}$$

A matrix $$A$$ is skew-symmetric if and only if $$A = \begin{pmatrix} \, 0 & z & -y \\ -z & 0 & x \\ y & -x & 0 \, \end{pmatrix}$$ Thus, for the sake of clarity, I introduce two functions: $$\overrightarrow{V} = \text{vector}(V)$$ which turns a skew-symmetric matrix $$V = \begin{pmatrix} \, 0 & V_3 & -V_2 \\ -V_3 & 0 & V_1 \\ V_2 & -V_1 & 0 \, \end{pmatrix}$$ into the vector column $$\overrightarrow{V} = \begin{pmatrix} V_1 \\ V_2 \\ V_3 \end{pmatrix}$$ and its inverse function $$V = \text{Matrix}(\overrightarrow{V})$$ which turns a vector column $$\overrightarrow{V}$$ into the corresponding skew-symmetric matrix $$V$$ with. We also will use the notation $$|V| = |\overrightarrow{V}| = \sqrt{V_1^2 + V_2^2 + V_3^2}$$ which is the length of the vector $$\overrightarrow{V}$$, and we can also think of it as the magnitude $$|V|$$ of the skew-symmetric matrix $$V$$.

The shape of the asteroid and the distribution of its mass is encoded in a $$3 \times 3$$ matrix $${\bf \mathcal{I}}$$, called the inertia tensor of the asteroid.

First, you have the coordinate system $$S = O\, \vec{e}_1\,\vec{e}_2\,\vec{e}_3$$ fixed in space, according to which you determine the positions of the point $$P$$ that fires upon the asteroid and the position of the center of mass of the asteroid $$A$$. Furthermore, you have another coordinate system $$S_A=A\, \vec{E}_1\,\vec{E}_2\,\vec{E}_3$$, firmly attached to the asteroid. Thus, the coordinate system $$S_A$$ rotates together with the asteroid. Despite the fact that the system $$S_A$$ is non-inertial, it is the right choice of system with respect to which to describe the dynamics of the motion because in $$S_A$$ the inertia tensor is constant, i.e. doesn't change with time, while in $$S$$ it does, making things complicated.

Thus, the idea is to study the (discrete) time evolution (i.e. evolution in small time-increments) of the quantities of motion in $$S_A$$ and then transform them into the system $$S$$. If $$\vec{X}$$ is the position vector of a point on the asteroid with respect to system $$S_A$$, the position vector $$\vec{x}(n)$$ at time $$n$$ of the same point with respect to $$S$$ can be calculate by applying the following transformation $$\vec{x}(n) = g(n) \,\vec{X} + {a}(n)$$ where $$n=0, 1, 2, 3, ...$$ is the discrete time, $$g(n) \in \text{SO}(3)$$ is a time-dependent orthogonal matrix, representing rotation at time $$n$$, and $$\vec{a}(n)$$ is the position of the center of mass $$A$$ of the asteroid at time $$n$$. Furthermore, let $$\vec{v}$$ be a fixed free vector with respect to $$S$$ (free vector is not a position vector). Then the same free vector, observed in $$S_A$$ will change with time $$n$$, so in $$S_A$$ vector $$\vec{v}$$ will be $$\vec{V}(n)$$ at time $$n$$ and the two are connected by $$\vec{v} = g(n) \, \vec{V}(n)$$ Next, after one time step (i.e. one time increment) time changes from $$n$$ to $$n+1$$ so we get $$\vec{v} = g(n+1)\, \vec{V}(n+1)$$ Thus, $$\vec{V}(n+1) = g(n+1)^{-1} \, \vec{v} = g(n+1)^{-1} \,g(n) \, \vec{V}(n)$$ Denote by $$W(n) = g(n)^{-1}\, g(n+1)\, \in \, \text{SO}(3)$$. Then $$\vec{V}(n+1) = W(n)^{-1} \, \vec{V}(n) = W(n)^T \, \vec{V}(n)$$ because $$W(n)^{-1} = W(n)^T$$ for any orthogonal matrix. The orthogonal operator $$W(n)^T$$ describes the one time step rotation of the vector $$\vec{V}(n)$$ to the vector $$\vec{V}(n+1)$$ in the coordinate system $$S_A$$. The matrix $$W(n)$$ is the discrete time analogue of angular velocity. If one finds the matrices $$W(n)$$ then one can generate the matrices $$g(n)$$ by $$g(n+1) = g(n) \, W(n)$$ If $$h>0$$ is the small time step (the time increment), then the actual time is $$t_{n} = n\, h$$ and $$t_{n+1} = t_{n} + h = (n+1)h$$. Then, during the time period $$h$$ between $$t_n$$ and $$t_{n+1}$$ the rotation $$W(n)$$ is around an axis determined by a vector $$\vec{\Omega}(n)$$ in $$S_A$$ for which $$\Omega(n) = \text{Matrix}(\vec{\Omega}(n))$$ and so \begin{align} W(n) = e^{h\, \Omega(n)} = {\bf 1} \, + \, \frac{\sin\Big(h\, \big|\Omega(n)\big|\Big)}{\big|\Omega(n)\big|}\,\, \Omega(n) \, + \, \frac{1 - \cos\Big(h\, \big|\Omega(n)\big|\Big)}{\big|\Omega(n)\big|^2}\,\, \Omega(n)^2\end{align} The vector $$\vec{\Omega}(n)$$ is the angular velocity written in the coordinate system $$S_A$$ and $${\bf 1}$$ is the identity matrix.

With analogy of the time continuous case of Euler's differential equations of motion, the latter are based on two laws of physics.

(i) The angular momentum $$\vec{m}$$ of the asteroid with respect to $$S$$ is constant, i.e. it doesn't change with time. This is the so called conservation of angular momentum law. Thus, $$\vec{m} = g(n) \, \vec{M}(n)$$ where $$\vec{M}(n)$$ is the same angular momentum but written with respect to the moving coordinate system $$S_A$$, so $$\vec{M}(n)$$ must change with time $$n$$. But for $$n+1$$ $$\vec{m} = g(n+1) \, \vec{M}(n+1)$$ Thus $$g(n) \, \vec{M}(n) = \vec{m} = g(n+1) \, \vec{M}(n+1)$$ which means that $$\vec{M}(n+1) = g(n+1)^{-1}\, g(n) \, \vec{M}(n)$$ and hence $$\vec{M}(n+1) = W(n)^T \, \vec{M}(n) = \Big( e^{-\,h\, \Omega(n)} \Big)\, \vec{M}(n)$$ where $$g(n)^{-1}\, g(n+1) = W(n) = e^{h\, \Omega(n)}$$ and $$W(n)^T = e^{-\, h\, \Omega(n)}$$.

(ii) The inertia matrix $$\mathcal{I}$$ gives the link between the angular velocity and the angular momentum $$\vec{M}(n) = \mathcal{I} \, \vec{\Omega}(n)$$ so the discrete time analogue of Euler's differential equations are the following discrete time difference equations

$$\vec{M}(n+1) = \left( e^{-\,h\, \text{Matrix}\big(\mathcal{I}^{-1} \, \vec{M}(n)\big)} \right)\, \vec{M}(n)$$ $$M(n) = \text{Matrix}(\vec{M}(n)) \,\,\, \text{ and } \,\,\, \vec{M}(n) = \text{Vector}({M}(n))$$ Generating $$\vec{M}(n)$$ leads to $$W(n) = e^{h\, \mathcal{I}^{-1} \, {M}(n)}$$ and consequently to $$g(n) = g_0\, W(1)\, W(2) \, ... \, W(n)$$

If $$Q$$ is a point on the asteroid with coordinates (i.e. position vector) $$\overrightarrow{Q} = \overrightarrow{AQ} = \begin{pmatrix}Q_1 \\ Q_2\\ Q_3 \end{pmatrix}$$ in the coordinate system $$S_A$$, the point $$Q$$ has coordinates (i.e. position vector) that change with time $$\vec{q}(n) = \overrightarrow{OQ}(n) = \begin{pmatrix} q_1(n) \\ q_2(n)\\ q_3(n) \end{pmatrix}$$ with respect to the coordinate system $$S$$ and $$\vec{q}(n) = g(n) \, \overrightarrow{Q} + \vec{a}(n) = \Big( g_0\, W(1)\, W(2) \, ... \, W(n) \Big) \, \overrightarrow{Q} + \vec{a}(n)$$

Algorithm

1. Start with initial position of asteroid's center of gravity $$A$$ given by a position vector $$\vec{a}(0)$$ with respect to the inertial coordinate system $$S = O\, \vec{e}_1\, \vec{e}_2\, \vec{e}_3$$.

2. Furthermore, fix an initial orientation of the asteroid with respect to $$S$$, given by a choice of three vectors forming the coordinate system $$S_A = A\, \vec{E}_1\, \vec{E}_2\, \vec{E}_3$$ attached to the asteroid at time $$n=0$$. The latter three vectors are written as $$3$$D vector columns with respect to $$S$$.

3. Form the orthogonal matrix $$g(0) = \Big(\vec{E}_1\, \vec{E}_2\, \vec{E}_3\Big)^T \, \in \, \text{SO}(3)$$.

4. Choose a point $$Q$$ on the asteroid, which is then given by a position vector $$\overrightarrow{Q}$$ with respect to $$S_A$$.

5. Take a diagonal inertia matrix $$\mathcal{I} = \text{diag}(I_1, I_2, I_3)$$.

6. Fix a time step $$h>0$$, say $$h = \frac{1}{25}$$ ($$25$$ frames per second is a reasonable choice to simulate smooth continuous motion on the screen, but the rate of change of frames can be chosen even faster if you want). Time is then $$t = n \, h$$ seconds, where $$n = 0, 1, 2, ...$$. Choose an integer $$N$$ so that time runs from $$t=0$$ to $$t = N\, h$$.

7. Give the asteroid an initial angular momentum $$\vec{M}(0)$$ written as a vector column with respect to $$S_A$$. Alternatively, you can choose initial angular velocity $$\vec{\Omega}(0) = \mathcal{I}^{-1} \, \vec{M}(0)$$.

8. Execute the following loop: For $$n = 0\, ... \, N$$ do the following \begin{align*} &\vec{\Omega}(n) = \mathcal{I}^{-1} \, \vec{M}(n) \\ &{\Omega}(n) = \text{Matrix}\Big( \vec{\Omega}(n) \Big) \\ &W(n) = e^{h\, {\Omega}(n)} = {\bf 1} \, + \, \frac{\sin\Big(h\, \big|\Omega(n)\big|\Big)}{\big|\Omega(n)\big|}\,\, \Omega(n) \, + \, \frac{1 - \cos\Big(h\, \big|\Omega(n)\big|\Big)}{\big|\Omega(n)\big|^2}\,\Omega(n)^2\\ &\vec{M}(n+1) = W(n)^T\, \vec{M}(n)\\ &g(n+1) = g(n) \, W(n) \end{align*}

The result is a sequence of matrices $$g(n) \, : \, n = 0\, ... \, N$$ and so the motion of point $$Q$$ with respect to the inertial non-moving frame $$S$$ is $$\vec{q}(n) = g(n) \, \overrightarrow{Q} + \vec{a}(n) = g(n) \, \overrightarrow{Q} + \vec{a}(0) + nh\, \vec{v}_A$$ where I have assumed that the center of mass $$A$$ of the asteroid moves uniformly with constant velocity $$\vec{v}_A$$ along a straight-line trajectory.

I will be very interested to get some feed back from you whether you have been able to realize the simulation of this kind of motion, regardless of the involvement of shooting projectiles at point $$Q$$ from the asteroid. So. please let me know whether it works.