From the IB Dynamics course, the equations of motion for a symmetric top can be written as:
$$ \frac{1}{2}I_1 \dot{\theta}^2 = E - \frac{(J_z-J_3 \cos \theta)^2}{2I_1 \sin^2 \theta} - mgh\cos\theta - \frac{J_3^2}{2I_3} $$$$ I_1\ddot{\theta} = mgh\sin\theta-\frac{J_3(J_z-J_3 \cos \theta)}{I_1 \sin \theta}+\frac{(J_z-J_3 \cos \theta)^2\cos\theta}{I_1 \sin^3\theta} $$$$ \dot{\phi} = \frac{J_z-J_3 \cos \theta}{I_1 \sin^2 \theta} $$$$ \dot{\chi} = \frac{J_3}{I_3}+\frac{J_3\cos^2\theta-J_z\cos\theta}{I_1 \sin^2\theta} $$By integrating the second equation, which is the time derivative of the first, to find $\dot{\theta}$, the angle $\theta(t)$ can be found and this can be used to integrate both $\chi$ and $\phi$.
First suitable values for the free parameters must be found. By writing the energy equation as the following:
$$ \frac{1}{2}I_1 \dot{\theta}^2 = E - \frac{(J_z-J_3 \cos \theta)^2}{2I_1 \sin^2 \theta} - mgh\cos\theta - \frac{J_3^2}{2I_3} = E - U_{eff}(\theta) $$The effective potential can be used to find the range of $\theta$ values for a given energy.
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.patches import FancyArrowPatch
from mpl_toolkits.mplot3d import proj3d
from matplotlib.animation import FuncAnimation
%matplotlib inline
from pycav import display
E = 71.
# I_1 and I_3
I = [0.1,0.05]
# J_z and J_3
J = [2.,2.5]
m = 1.
h = 1.
def Ueff(theta,I,J,m,h):
return (J[0]-J[1]*np.cos(theta))**2/(2*I[0]*(np.sin(theta))**2)+m*9.81*h*np.cos(theta)+J[1]**2/(2*I[1])
By looking at the plot below, one can find the maximum and minimal values of $\theta$ for the given energy value, $E$.
theta = np.linspace(0.2,np.pi/2.,1000)
fig1 = plt.figure(figsize = (15,9));
ax1 = fig1.add_subplot(111);
ax1.plot(theta,Ueff(theta,I,J,m,h));
ax1.plot(theta,(J[0]-J[1]*np.cos(theta))**2/(2*I[0]*(np.sin(theta))**2)+J[1]**2/(2*I[1]));
ax1.plot(theta,m*9.81*h*np.cos(theta));
ax1.axhline(E,c = 'k',ls = '-.');
ax1.set_xlim((theta.min(),theta.max()));
ax1.set_xlabel('Theta');
ax1.legend(('Effective Potential','Ang. Mom. Term','Gravitational Term','Energy'));
With suitable free parameters chosen, the equations can be integrated using SciPy's odeint function to give us the angles as functions of time
Here we define $\omega = \dot{\theta}$ so the set of first order differential equations we are solving are:
$$ \dot{\omega} = \frac{1}{I_1} \left( mgh\sin\theta-\frac{J_3(J_z-J_3 \cos \theta)}{I_1 \sin \theta}+\frac{(J_z-J_3 \cos \theta)^2\cos\theta}{I_1 \sin^3\theta} \right) $$$$ \dot{\theta} = \omega $$$$ \dot{\phi} = \frac{J_z-J_3 \cos \theta}{I_1 \sin^2 \theta} $$$$ \dot{\chi} = \frac{J_3}{I_3}+\frac{J_3\cos^2\theta-J_z\cos\theta}{I_1 \sin^2\theta} $$odeint takes a list of starting values, [$\omega_0$, $\theta_0$, $\phi_0$, $\chi_0$], and takes a function which returns a list of the time derivatives, [$\dot{\omega}$, $\dot{\theta}$, $\dot{\phi}$, $\dot{\chi}$], and returns a list containing the input variables evaluated at the times given, [$\omega(t)$, $\theta(t)$, $\phi(t)$, $\chi(t)$]
def omega_dot(theta,E,I,J,m,h):
return (1/I[0])*(m*9.81*h*np.sin(theta)-
(J[1]*(J[0]-J[1]*np.cos(theta)))/(I[0]*np.sin(theta))+
(J[0]-J[1]*np.cos(theta))**2*np.cos(theta)/(I[0]*(np.sin(theta))**3))
def theta_dot(theta,E,I,J,m,h):
rot_energy = E-Ueff(theta,I,J,m,h)
return np.sqrt(2*rot_energy/I[0])
def phi_dot(theta,E,I,J,m,h):
return (J[0]-J[1]*np.cos(theta))/(I[0]*(np.sin(theta))**2)
def chi_dot(theta,E,I,J,m,h):
return (J[1]*(np.cos(theta))**2-J[0]*np.cos(theta))/(I[0]*(np.sin(theta))**2)+J[1]/I[1]
def derivatives(q,t,E,I,J,m,h):
omega = q[0]
theta = q[1]
dqdt = [omega_dot(theta,E,I,J,m,h),omega,phi_dot(theta,E,I,J,m,h),chi_dot(theta,E,I,J,m,h)]
return dqdt
starting_angles = [0.8,0.,0.]
q0 = [theta_dot(starting_angles[0],E,I,J,m,h),starting_angles[0],starting_angles[1],starting_angles[2]]
t = np.linspace(0.,3.5,500)
sol = odeint(derivatives,q0,t,args = (E,I,J,m,h))
Converting these angles into Cartesian coordinates allows us to plot the trajectory of the centre of the mass
theta = sol[:,1]
phi = sol[:,2]
chi = sol[:,3]
x,y,z = h*np.sin(theta)*np.cos(phi),h*np.sin(theta)*np.sin(phi),h*np.cos(theta)
# Create cubic bounding box to simulate equal aspect ratio
max_range = np.array([x.max()-x.min(), y.max()-y.min(), z.max()-z.min()]).max()
Xb = 0.5*max_range*np.array((-1.,1.)) + 0.5*(x.max()+x.min())
Yb = 0.5*max_range*np.array((-1.,1.)) + 0.5*(y.max()+y.min())
Zb = 0.5*max_range*np.array((-1.,1.)) + 0.5*(z.max()+z.min())
fig2 = plt.figure(figsize = (18,8));
fig2.suptitle('Resultant Motion of Gyroscope', fontsize = 15);
ax2 = fig2.add_subplot(121, projection='3d');
ax2.set_aspect('equal');
ax2.set_xlim3d((Xb[0],Xb[-1]));
ax2.set_ylim3d((Yb[0],Yb[-1]));
ax2.set_zlim3d((Zb[0],Zb[-1]));
ax2.plot(x,y,z,'r');
ax3 = fig2.add_subplot(122);
ax3.plot(t,x,'r');
ax3.plot(t,y,'b');
ax3.plot(t,z,'k');
ax3.legend(('x','y','z'));
ax3.set_xlabel('t');
class Arrow3D(FancyArrowPatch):
def __init__(self, xs, ys, zs, *args, **kwargs):
FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
self._verts3d = xs, ys, zs
def draw(self, renderer):
xs3d, ys3d, zs3d = self._verts3d
xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
FancyArrowPatch.draw(self, renderer)
fig = plt.figure(figsize = (12,12));
ax = fig.add_subplot(111, projection='3d');
ax.set_aspect('equal');
qx = 0.5*x+0.1*h*(np.cos(theta)*np.cos(phi)*np.cos(chi)-np.sin(phi)*np.sin(chi))
qy = 0.5*y+0.1*h*(np.cos(theta)*np.sin(phi)*np.cos(chi)+np.cos(phi)*np.sin(chi))
qz = 0.5*z+0.1*h*(-np.sin(theta)*np.cos(chi))
def nextframe(arg):
n = arg*2
ax.clear()
ax.plot([0.5*x[n],qx[n]],[0.5*y[n],qy[n]],[0.5*z[n],qz[n]],'k-',lw = 3)
a = Arrow3D([0., x[n]], [0., y[n]],
[0., z[n]], mutation_scale=20,
lw=3, arrowstyle="-|>", color="r")
ax.add_artist(a)
ax.plot(x[:n],y[:n],z[:n],'r')
ax.set_xlim3d((Xb[0],Xb[-1]))
ax.set_ylim3d((Yb[0],Yb[-1]))
ax.set_zlim3d((Zb[0],Zb[-1]))
animate = FuncAnimation(fig,nextframe, interval = 75, frames = int(len(t)/2), repeat = False)
animate = display.create_animation(animate, temp = True)
display.display_animation(animate)