Solving a Coupled System of First Order, Nonlinear Differential Equations to Model Rotation on a Changing Axis I am looking to model the rotation of some vector representing a rigid body about its origin given:


*

*the initial vector $\vec x_0$

*a constant $\phi$ denoting the magnitude of the initial angular velocity in rads/s

*a constant $\lambda$ denoting the magnitude of the constant angular acceleration in rads/s^2

*a constant unit vector $\hat v_0$ denoting the initial axis of rotation on with BOTH the angular velocity and acceleration are applied; its components are denoted by: $\hat v_x \hat v_y \hat v_z$
After time $t_0$, the new axis of rotation is the vector $\hat v_0$ rotated into the new orientation of the rigid body at time t. Basically, if the rotation of the body at any time t could (even though it clearly cannot) be modeled by some rotation matrix M(t), then the instantaneous axis vector would be $\hat v_0$ transformed by M(t).
Searching for a solution to this problem, I came across the following paper: https://arxiv.org/pdf/1604.08139v1.pdf
In section 3, it outlines a method of modeling a rotation of the type I see to model using quaternions. The rotated version of $\hat x_0$ would be given by
\begin{equation}
  \mathbf{x}'(t) = R(t)\, \mathbf{x}\, \bar{R}(t),
\end{equation}
where R(t) is some quaternion used to apply the accumulated rotation from time 0 to time t, and $\mathbf{x}$ is a quaternion with a zero scaler part and vector part of $\vec x_0$.
The paper arrives at the following differential equation to model R(t) in terms of w(t), the angular velocity vector with respect to time(direction is axis, magnitude is magnitude):
\begin{equation}
  \dot{R}(t) = \frac{1}{2}\, \boldsymbol{\omega}(t)\, R(t).
\end{equation}
This models the case where w(t) is given in terms of the "global axes". It is not relative to the current orientation of the rigid body. However, to suite my purposes, w(t) must be measured in respect to the body. w(t) in terms of the angular velocity with respect to the body is given by:
$\boldsymbol{\omega}(t) = R(t)\, \boldsymbol{\omega}'(t)\, \bar{R}(t)$
Substituting, the aforementioned differential equation becomes:
\begin{equation}
  \dot{R}(t) = \frac{1}{2}\, R(t)\, \boldsymbol{\omega}'(t).
\end{equation}
The quaternion R(t) may by decomposed into four parametric functions:
[$q_w(t)$, $q_x(t)$, $q_y(t)$, $q_z(t)$]
Based on the initial givens, $\boldsymbol{\omega}'(t)$ may also be written as follows: [0, ($\lambda$t+$\phi$)$\hat v_0$]
Applying rules of quaternion multiplication on the final differential equation, I arrive at the following system of coupled, first-order differential equations:
\begin{align}
\dot{q}_w(t) &= -1/2\,(\lambda\, t+\phi)(q_x(t)\hat v_x + q_y(t)\hat v_y + q_z(t)\hat v_z), \\
\dot{q}_x(t) &= 1/2\,(\lambda\, t+\phi)(q_w(t)\hat v_x + q_y(t)\hat v_z - q_z(t)\hat v_y), \\
\dot{q}_y(t) &= 1/2\,(\lambda\, t+\phi)(q_w(t)\hat v_y + q_z(t)\hat v_x - q_x(t)\hat v_z), \\
\dot{q}_z(t) &= 1/2\,(\lambda\, t+\phi)(q_w(t)\hat v_z + q_x(t)\hat v_y - q_y(t)\hat v_x).
\end{align}
Is this system solvable analytically? If not, would removing angular acceleration ($\lambda$) make any difference? 
This is not a homework question or anything. I'm just pondering it and was hoping someone could help or lend insight. Any assistance would be greatly appreciated.  
 A: If the angular velocity is just some scalar function of the time times a constant vector you are essentially only rotating around that vector. The angle rotated around that vector would just be the integral of that scalar function. For ease of notation I will denote a quaternion as $(q_w,\vec{q}_v)$ with $\vec{q}_v=(q_x,q_y,q_z)$. In your case with an initial condition of the quaternion of $R(0)=(1,\vec{0})$ the analytical solution would be
$$
R(t) = \left(
\cos\left(\frac{\lambda\,t^2 + 2\,\phi\,t}{4}\right),
\hat{v}_0 \sin\left(\frac{\lambda\,t^2 + 2\,\phi\,t}{4}\right)
\right). \tag{1}
$$
In order to show that this is an analytical solution one can differentiate $(1)$ and show that that is equivalent to plugging $(1)$ into the differential equation. Differentiating $(1)$ gives
$$
\dot{R}(t) = \frac{\lambda\,t + \phi}{2} \left(
-\sin\left(\frac{\lambda\,t^2 + 2\,\phi\,t}{4}\right),
\hat{v}_0 \cos\left(\frac{\lambda\,t^2 + 2\,\phi\,t}{4}\right)
\right). \tag{2}
$$
Plugging $(1)$ into the differential equation gives
$$
\dot{R}(t) = \frac{\lambda\,t + \phi}{2} \left(
-\hat{v}_0^\top \hat{v}_0 \sin\left(\frac{\lambda\,t^2 + 2\,\phi\,t}{4}\right),
\hat{v}_0 \cos\left(\frac{\lambda\,t^2 + 2\,\phi\,t}{4}\right) + 
\hat{v}_0 \sin\left(\frac{\lambda\,t^2 + 2\,\phi\,t}{4}\right) \times \hat{v}_0
\right). \tag{3}
$$
When using that $\hat{v}_0^\top \hat{v}_0 = 1$ and $\hat{v}_0 \times \hat{v}_0=\vec{0}$ it can be shown that $(2)$ is equivalent to $(3)$.
If $R(0)\neq(1,\vec{0})$ one could use $R(t)=R(0)\,R'(t)$ with $R'(t)$ defined using $(1)$. Namely by definition $R'(0)=(1,\vec{0})$ and $\dot{R}'(t) = \bar{R}(0)\,\dot{R}(t) = 1/2\,\bar{R}(0)\,R(t)\,\omega'(t) = 1/2\,R'(t)\,\omega'(t)$.
