SHARE
TWEET

Untitled

a guest Feb 23rd, 2019 62 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. g = 9.81
  2. l = 1.6
  3. l_big = 2.0
  4. l_small = 1.6
  5. m = 0.01
  6. alpha = 0.4
  7. k = 100
  8. def sh2(r1,t):
  9.  
  10.     theta1,omega1 = r1
  11.     sh2_theta1 = omega1/(l + ((1/2)*alpha*(1+np.tanh(theta1*omega1*k))))**2
  12.  
  13.     sh2_omega1 = -g*(l + ((1/2)*alpha*(1+np.tanh(theta1*omega1*k))))*sin(theta1)
  14.  
  15.     #return sh2_theta1, sh2_omega1
  16.  
  17.     return np.array([sh2_theta1, sh2_omega1],float)
  18.  
  19. init_state = np.radians([69.0,0])
  20. dt = 1/40
  21. time = np.arange(0,10.0,dt)
  22. timexo = np.arange(0,10.0,dt)
  23.  
  24. state2 = integrate.odeint(sh2,init_state,time)
  25. #print(state2)
  26. print(len(state2),len(timexo))
  27.  
  28. plt.plot(timexo[0:2500],state2[0:2500])
  29. plt.show()
  30.  
  31. #phase plot attempt
  32. # initial values
  33. x_0 = 0 # intial angular position
  34. v_0 = 1.0 # initial angular momentum
  35. t_0 = 0 # initial time
  36.  
  37. # initial y-vector from initial position and momentum
  38. y0 = np.array([x_0,v_0])
  39.  
  40.  
  41. # max value of time and points in time to integrate to
  42. t_max = 10
  43. N_spacing_in_t = 1000
  44.  
  45. # create vector of time points you want to evaluate
  46. t = np.linspace(t_0,t_max,N_spacing_in_t)
  47.  
  48. # create vector of positions for those times
  49. y_result = np.zeros((len(t), 2))
  50.  
  51. # loop through all demanded time points
  52. for it, t_ in enumerate(t):
  53.  
  54.     # get result of ODE integration up to the demanded time
  55.     y = integrate.odeint(sh2,init_state,t_)
  56.  
  57.     # write result to result vector
  58.     y_result[it,:] = y
  59.  
  60. # get angle and angular momentum
  61. angle = y_result[:,0]
  62. angular_momentum = y_result[:,1]
  63.  
  64. # plot result
  65. pl.plot(angle, angular_momentum,'-',lw=1)
  66. pl.xlabel('angle $x$')
  67. pl.ylabel('angular momentum $v$')
  68.  
  69. pl.gcf().savefig('pendulum_single_run.png',dpi=300)
  70.  
  71. pl.show()
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top