In [9]:
from __future__ import division, print_function
from ivisual import *
from numpy import *
import matplotlib.pyplot as plt
In [13]:
k=5 #stiffness
m=0.05 #mass
b=0.022
omega0=sqrt(k/m) #natural frequency
beta=b/2/m #decay of amplitude
omega1=sqrt(omega0*omega0-beta*beta) #oscillation frequency
delta=-atan(-beta/omega1) #phase
A0=0.1 #constant
x=A0*cos(-delta) #initial position; technically delta and A0 depend on this value
vx=0 #initial velocity; this was used to calculate delta
t=0
dt=0.0001
tlist=[]
xlist=[]
vxlist=[]
extremalist=[]
while t<4:
vxi=vx
Fnet=-k*x-b*vx
vx=vx+Fnet/m*dt
x=x+vx*dt
t=t+dt
tlist.append(t)
xlist.append(x)
vxlist.append(vx)
if((vxi>0 and vx<0) or (vxi<0 and vx>0) ):
extremalist.append(t)
print("extremum at t=",t)
In [14]:
%matplotlib inline
In [15]:
#analytic function
xtheor=exp(-beta*array(tlist))*A0*cos(omega1*array(tlist)-delta)
plt.title('x vs. t')
plt.xlabel('t (s)')
plt.ylabel('x (m)')
plt.plot(tlist,xlist,'m.', tlist,xtheor,'y.')
plt.show()
In [17]:
#theoretical extrema
for n in arange(1,13,1):
print("extremum at t=",(n*pi)/omega1)
In [ ]: