Note: These tests are primarily to ensure that I understand what is going on, and things are working reasonably. I have checked that the results converge appropriately. However, this is not the preferred method for finding a frame from an angular velocity. Because the data are fixed at specified intervals, the accuracy is not great. If possible, it would be better to integrate the system as an ODE. Such a procedure is demonstrated in the C++ code TestAngularVelocityIntegration.cpp in this directory.
Because of the discrete nature of the data used to do the integration -- and in particular, because interpolation and differentiation near the beginning and end of such data is not very accurate -- the beginning and ends of the integrated frame data are not very accurate. We remove these errors by forcing the analytical solution and the numerical solution to be identical ~20 steps into the data. If they are the same later, we have a good solution.
This test corresponds to a rotation frequency of roughly $\alpha$ from $T_i$ to $T_f$ with time steps of $dt$ (each set below).
In [2]:
alpha = 0.05
T_i = 0.0
T_f = 2000.0
dT = 0.1
cushion = 20 # Number of indices to avoid at the beginning and end
import Quaternions
# Construct a frame R_1
T = arange(T_i, T_f, dT)
R_1 = array([ Quaternions.exp(Quaternions.xHat*(alpha*t/2.))
*Quaternions.exp(Quaternions.yHat*((alpha*t**1.2)/2.))
*Quaternions.exp(Quaternions.zHat*((alpha*t**1.3)/2.))
for t in T])
# Get the angular velocity Omega of that frame
Rdot_1 = array([ Quaternions.xHat*R_1[i]*(alpha/2.)
+(Quaternions.exp(Quaternions.xHat*(alpha*t/2.))
*(1.2*alpha*t**0.2/2.))*Quaternions.yHat
*Quaternions.exp(Quaternions.yHat*((alpha*t**1.2)/2.))
*Quaternions.exp(Quaternions.zHat*((alpha*t**1.3)/2.))
+R_1[i]*Quaternions.zHat*((1.3*alpha*t**0.3)/2)
for i,t in enumerate(T)])
Omega = array([Rdot*Quaternions.conjugate(R)*2 for Rdot,R in zip(Rdot_1, R_1)])
# Integrate a new frame R_2 from Omega
%time R_2 = array(Quaternions.FrameFromAngularVelocity(Omega, T))
# Check to make sure that R_1 and R_2 are the same
R_Delta = Quaternions.RDelta(R_1, R_2, cushion)
theta_Delta = array(Quaternions.angle(R_Delta))
semilogy(T, theta_Delta);
print(theta_Delta[-cushion])
# Error goes approximately as t**2
def f(t):
return t**2
semilogy(T, theta_Delta[-cushion]*f(T)/f(T[-cushion]));
For certain cases, we want to restrict the quaternion logarithm to a 2-D subspace. For example, if we have a unit vector $\hat{v}$ and have an angular velocity vector $\vec{\Omega}$, we may want to construct a rotor $R$ such that its frame has angular velocity $\vec{\Omega}$, and it gives the unit vector as $\hat{v} = R\, \hat{z}\, \bar{R}$, for example. This is almost the same situation as above, but not quite.
If all we care about is the direction of that one vector, we only need two degrees of freedom, so integrating the three degrees of freedom in the logarithm of $R$ is too much. Instead, we can restrict $\log R$ to a 2-D subspace, and integrate to find it.
Here again, the solution will be inaccurate near the beginning and end of the integrated data because the errors interpolating and differentiating discrete data are largest near the beginning and end. But this time, it's significantly harder to get rid of these errors. Above, we relied on the fact that the angular velocities of the two frames had to be the same, so there was only one way to adjust the numerical frame while preserving that condition. Here, there is no such simplification. Instead, we simply subtract the magnitude of the error after several steps from the magnitude of the error at other times. There may be additional error in the direction of the vector.
In [3]:
alpha = 0.05
T_i = 0.0
T_f = 2000.0
dT = 0.1
import Quaternions
# Construct a frame R_1
T = arange(T_i, T_f, dT)
R_1 = array([ Quaternions.exp(Quaternions.xHat*(alpha*t/2.))
*Quaternions.exp(Quaternions.yHat*((alpha*t**1.2)/2.))
*Quaternions.exp(Quaternions.zHat*((alpha*t**1.3)/2.))
for t in T])
# Get the angular velocity Omega of that frame
Rdot_1 = array([ Quaternions.xHat*R_1[i]*(alpha/2.)
+(Quaternions.exp(Quaternions.xHat*(alpha*t/2.))
*(1.2*alpha*t**0.2/2.))*Quaternions.yHat
*Quaternions.exp(Quaternions.yHat*((alpha*t**1.2)/2.))
*Quaternions.exp(Quaternions.zHat*((alpha*t**1.3)/2.))
+R_1[i]*Quaternions.zHat*((1.3*alpha*t**0.3)/2.)
for i,t in enumerate(T)])
Omega = array([Rdot*Quaternions.conjugate(R)*2 for Rdot,R in zip(Rdot_1, R_1)])
# Integrate a new frame R_2 from Omega
%time R_2 = array(Quaternions.FrameFromAngularVelocity_2D(Omega, T))
#R_2 = R_2 * (R_2[cushion].inverse()*R_1[cushion])
#R_2 = R_2 * Quaternions.exp(Quaternions.zHat*())
# Integrate a new frame R_3 from Omega
%time R_3 = array(Quaternions.FrameFromAngularVelocity(Omega, T))
R_3 = R_3 * (R_3[cushion].inverse()*R_1[cushion])
# Test whether or not R_1 and R_2 give rise to the same z vector
z_1 = array([R*Quaternions.zHat*R.conjugate() for R in R_1])
z_2 = array([R*Quaternions.zHat*R.conjugate() for R in R_2])
z_3 = array([R*Quaternions.zHat*R.conjugate() for R in R_3])
z_Delta12 = z_1 - z_2
zDelta12Mag = array([abs(d.abs()-z_Delta12[cushion].abs()) for d in z_Delta12])
semilogy(T, zDelta12Mag, label='Using 2-D integration');
z_Delta13 = z_1 - z_3
zDelta13Mag = array([d.abs() for d in z_Delta13])
semilogy(T, zDelta13Mag, label='Using 3-D integration');
#xlim((0, 2*cushion));
legend(loc='upper left');
print(zDelta12Mag[-cushion],zDelta13Mag[-cushion])
Not only is the 2-D integration more accurate and stable, it is also significantly faster than the 3-D integration.