How to Use Quaternions for Integrals, Differential Equations, Physics and Robotics?

Question

For the case of translation, it is very simple to use Euler's method (among other methods that are far more efficient) to interpolate between positions with increments for time dt while applying forces, torques, and everything. You can keep track of Energies and momenta, integrate things, use numerical methods, and whatnot. Everything is "simple".

I am trying to model the flapping of birds with 3 degrees of rotational freedom. As it turns out, rotation is non-commutative. So, it isn't really that simple. Using simple shifts in coordinates to Euler angles and treating them like x,y, and z doesn't really work (or maybe I just haven't figured out how to make it work without running into huge errors).

SLERP smoothing is what I came across which uses quaternions to more conveniently do such math with rotating rigid bodies.

The way I can just take Work to be the integral of F dot dx for all coordinates x,y, and z; how do I use quaternions to do such integrals?

https://eaglepubs.erau.edu/introductiontoaerospaceflightvehicles/chapter/finite-wing-characteristics/ How, for instance, do integrals work when rotations are represented as quaternions?

Say I rotate from a state (r1,p1,y1) to (r2,p2,y2). What would the lift force (or any quantity that I should calculate as a function of time, dx, dy, dz, r, theta, or any other quantity) be for a certain time interval? Or something like that - any integral.

From what I understood, for rotations, you don't add a difference vector; instead, you multiply ("apply") a rotational transformation. How do differential and integral calculus work, then? How do I solve differential equations? How do I convert from one frame to another while the motion is in one frame, but I only have equations of properties in the other?

Answer 1

To do physics/calculus in the rotational domain you generally want to follow the procedure below. Note that this assumes a free body, with no external constraints (i.e. joints to other bodies) as that is a whole other topic. Not to say that is doesn't apply to your simulation problem but it is not the topic of the OP.

1. Calculate torque on the body $\tau$ in the local frame.
2. Solve for rotational acceleration in the local frame $\alpha = \mathbf I^{-1}\tau$ where $\mathbf I$ is the rotational inertia of the body expressed in the local frame.
3. Because $\alpha$ is a vector, you can integrate it in what ever way you want to (i.e. Runge-Kutta) to get the angular velocity $\omega$ in the local frame.
4. key point calculate the time derivative of the quaternions from the angular velocity using $\dot{q} =\frac{1}{2}q\omega$. Where $\omega$ is expressed as a pure quaternion (i.e. the real part is 0, and the imaginary part is the angular velocity.
5. Integrate $\dot{q}$ in whatever way you want to to get the new $q$.
6. Normalize $q$ to a magniqude of 1.
7. Convert $q$ to what every representation you like (i.e. rotation matrix, Euler angles, etc.)

here are some good resources:

https://ahrs.readthedocs.io/en/latest/filters/angular.html

https://www.euclideanspace.com/physics/kinematics/angularvelocity/

Answer 2

Suppose $\hat{M}$ is the matrix of multiplications to apply to state vector $v$ to update for a change of time, $\Delta t$. We usually find/expect that for small $\Delta t$, $\hat{M}$ is close to the identity matrix, $I$, (that is, we only change the coordinates of $v$ in whatever coordinate system by a small amount). Then it is convenient to define $M$ by $$ M = \hat{M} - I \text{.} $$ Then, $$ v(t+\Delta t) = (I+M)v(t) \text{.} $$

Notice that we can make a "better" step by taking a half step twice. When $M$ is "small", \begin{align*} (I+M/2)^2 &= I^2 + 2 M/2 + M^2/4 \\ &= I + M + M^2/4 \end{align*} If splitting in half is good, what happens if we split into $n$ (sub-)steps? (We apply Newton's binomial formula to obtain...) \begin{align*} (I+M/n)^n &= I^n + n I^{n-1}M/n + \frac{n(n-1)}{2} I^{n-2}(M/n)^2 + \cdots \\ &\qquad + n I (M/n)^{n-1} + (M/n)^n \\ &= I + M + \frac{n-1}{2n} M^2 + \cdots + M^{n-1}/n^{n-2} + (M/n)^n \text{,} \end{align*} What if we split into infinitely many sub-steps? $$ \lim_{n \rightarrow \infty} (I+M/n)^n = \mathrm{e}^{M} \text{.} $$ (As seen at https://math.stackexchange.com/a/722309/123905 )

So the "best" update for $v$ is $$ v(t+\Delta t) = \mathrm{e}^M v(t) \text{.} $$ (Conveniently, we can also write $$ v\left(t + \frac{1}{n} \Delta t\right) = \mathrm{e}^{M/n} v(t) $$ if we should need to zoom in the time steps to take smaller steps (to match the time steps of some other process that we are suddenly interacting with, for instance).)

Note, however, that we have written an integral in disguise. Apply $M$ for $\Delta t$ time, starting at $v$ and you'll obtain $v(t+\Delta t)$ -- that is, the integral of the effect of $M$ applied for duration $\Delta t$, ..., and we might recognize the differential equation that gives this solution. $$ y' = k y $$ has solution $$ y = \mathrm{e}^{kt} C \text{,}$$ where $C$ is the constant of integration (that selects which of the family of solutions happens to pass through our initial condition, the role played by $v$ in our matrix equation). That is, translating back to our matrix notation $$ v(t+\Delta t) = \mathrm{e}^{(M/\Delta t)\Delta t} v(t) \text{ and } $$ $$ v'(t) = (M/\Delta t) v(t) \text{.} $$

Feel free to ask follow-up questions, or modify the Question to more precisely highlight what is missing here (and then, maybe, let we answerers know that the Question has changed by commenting on the Answers).