[SciPy-user] Bicycle Wheel Precession

[Ryan Krauss]

I am trying to teach my dynamics students about teh precession of a
bicycle wheel like in this video:

http://commons.bcit.ca/physics/video/precession.shtml

If the x axis is along the axle of the wheel and the z axis points up,
I think the motion can be modeled using the following ode:

def NewtonEuler(x,t,M,I):
    Ixx=I[0]
    Iyy=I[1]
    Izz=I[2]
    wx=x[0]
    wy=x[1]
    wz=x[2]
    dxdt=zeros(3,)
    dxdt[0]=((Iyy-Izz)*wy*wz+M[0])/Ixx
    dxdt[1]=((Izz-Ixx)*wx*wz+M[1])/Iyy
    dxdt[2]=((Ixx-Iyy)*wx*wy+M[2])/Izz
    return dxdt

I use this function along with integrate.odeint in the attached file,
but the surprising thing is that wz is sinusoidal rather than nearly
constant.  Experimentally, I see nearly constant precession.

Can anyone help me find what is wrong with my approach?

Thanks,

Ryan
-------------- next part --------------
from scipy import *
from pylab import figure, cla, clf, plot, subplot, show, ylabel, xlabel, xlim, ylim, semilogx, legend, title, savefig, yticks, grid, rcParams

#from  IPython.Debugger import Pdb

import copy, os, sys


R=0.5#meters
m=5.0;#kg

Ixx=0.5*m*R**2
Iyy=0.25*m*R**2
Izz=Iyy

wx=2.0*pi#rad/sec roughly

def NewtonEuler(x,t,M,I):
    Ixx=I[0]
    Iyy=I[1]
    Izz=I[2]
    wx=x[0]
    wy=x[1]
    wz=x[2]
    dxdt=zeros(3,)
    dxdt[0]=((Iyy-Izz)*wy*wz+M[0])/Ixx
    dxdt[1]=((Izz-Ixx)*wx*wz+M[1])/Iyy
    dxdt[2]=((Ixx-Iyy)*wx*wy+M[2])/Izz
    return dxdt


def omegaplot(fi, t, out, mylegend=['$\\omega_x$','$\\omega_y$','$\\omega_z$']):
    figure(fi)
    cla()
    for col in out.T:
        plot(t,col)
    xlabel('Time (sec)')
    ylabel('Angular Velocity (rad/sec)')
    if mylegend:
        legend(mylegend)


wo=[wx,0.0,0.0]

I=[Ixx,Iyy,Izz]

d=0.05
g=9.81
My=d*g*m
M=[0,My,0]

t=arange(0,5,0.001)

out=integrate.odeint(NewtonEuler, wo, t, args=(M,I,))

figure(1)
cla()
plot(t,out[:,0],t,out[:,1],t,out[:,2])
legend(['$\\omega_x$','$\\omega_y$','$\\omega_z$'])
xlabel('Time (sec)')
ylabel('Angular Velocity (rad/sec)')
savefig('wx_positive.eps')

############
w2=[-wx,0.0,0.0]
out2=integrate.odeint(NewtonEuler, w2, t, args=(M,I,))

figure(2)
cla()
plot(t,out2[:,0],t,out2[:,1],t,out2[:,2])
legend(['$\\omega_x$','$\\omega_y$','$\\omega_z$'])
xlabel('Time (sec)')
ylabel('Angular Velocity (rad/sec)')
savefig('wx_negative.eps')


## def NewtonEulerSO(x,t,M,I):
##     Ixx=I[0]
##     Iyy=I[1]
##     Izz=I[2]
##     thx=x[0]
##     thy=x[1]
##     thz=x[2]
##     wx=x[3]
##     wy=x[4]
##     wz=x[5]
##     dxdt=zeros(6,)
##     dxdt[0]=wx
##     dxdt[1]=wy
##     dxdt[2]=wz
##     dxdt[3]=((Iyy-Izz)*wy*wz+M[0])/Ixx
##     dxdt[4]=((Izz-Ixx)*wx*wz+M[1]*cos(thx))/Iyy
##     dxdt[5]=((Ixx-Iyy)*wx*wy+M[2])/Izz
##     return dxdt


## Xo=[0.0,0.0,0.0,wx,0.0,0.0]

## outSO=integrate.odeint(NewtonEulerSO, Xo, t, args=(M,I,))

## xyzlist=['x','y','z']
## omegaleg=['$\\omega_'+item+'$' for item in xyzlist]
## thetalist=['$\\theta_'+item+'$' for item in xyzlist]
## solegend=thetalist+omegaleg
## omegaplot(3,t,outSO,mylegend=solegend)

show()