Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #No.11 #Calculate the period of simple harmonic oscillator
- import numpy as np
- from pylab import plot, show, grid, xlabel, ylabel
- import matplotlib.pyplot as plt
- def main ():
- #Total time
- T = 10.0 #how far the graph will go in seconds
- #number of steps
- N = 100
- #Time step size in sec
- E = 0.1
- #create array with steps+time step value from 0to time final
- #t = np.linspace(0.0,t,N+1)
- #print t
- #allocating arrays for velocity and position
- v=np.zeros(N+1)
- x=np.zeros(N+1)
- a=np.zeros(N+1)
- t=np.linspace(0.0, T, N+1)
- rev=0
- period=0
- #set the initial value for velocity and position
- v[0]=0.0
- x[0]=1.0
- a[0]=-1.0
- v[1]=v[0]+a[0]*E/2
- t[0]=0.0
- #looping for position and velocity
- for i in range(1,N):
- t[i]=t[i-1]+E
- x[i]=x[i-1]+E*v[i]
- a[i]=-x[i]
- v[i+1]=v[i]+E*a[i]
- if x[i]*x[i-1] < 0:
- if rev == 0:
- time_old=t[i]
- else:
- period[rev]=2*(t(i)-time_old)
- print period[rev]
- time_old=t[i]
- #rev=rev+1
- #print period[rev]
- period_mean=np.mean(period[rev])
- error=np.std(period[rev])/np.sqrt(rev)
- print ('Average period = ',period_mean, error)
- #plotting graph position vs time
- plt.figure(1)
- plt.plot (t, x, 'r')
- plt.xlabel('time', fontsize=16)
- plt.ylabel('position', fontsize=16)
- plt.grid(True)
- plt.show()
- if __name__ == "__main__":
- main()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement