双摆问题

这个动画说明了双摆问题。

双摆公式从 http://www.physics.usyd.edu.au/~wheat/dpend_html/solve_dpend.c 的C代码翻译而来。 

  1. from numpy import sin, cos
  2. import numpy as np
  3. import matplotlib.pyplot as plt
  4. import scipy.integrate as integrate
  5. import matplotlib.animation as animation
  7. G = 9.8  # acceleration due to gravity, in m/s^2
  8. L1 = 1.0  # length of pendulum 1 in m
  9. L2 = 1.0  # length of pendulum 2 in m
  10. M1 = 1.0  # mass of pendulum 1 in kg
  11. M2 = 1.0  # mass of pendulum 2 in kg
  14. def derivs(state, t):
  16.     dydx = np.zeros_like(state)
  17.     dydx[0] = state[1]
  19.     del_ = state[2] - state[0]
  20.     den1 = (M1 + M2)*L1 - M2*L1*cos(del_)*cos(del_)
  21.     dydx[1] = (M2*L1*state[1]*state[1]*sin(del_)*cos(del_) +
  22.                M2*G*sin(state[2])*cos(del_) +
  23.                M2*L2*state[3]*state[3]*sin(del_) -
  24.                (M1 + M2)*G*sin(state[0]))/den1
  26.     dydx[2] = state[3]
  28.     den2 = (L2/L1)*den1
  29.     dydx[3] = (-M2*L2*state[3]*state[3]*sin(del_)*cos(del_) +
  30.                (M1 + M2)*G*sin(state[0])*cos(del_) -
  31.                (M1 + M2)*L1*state[1]*state[1]*sin(del_) -
  32.                (M1 + M2)*G*sin(state[2]))/den2
  34.     return dydx
  36. # create a time array from 0..100 sampled at 0.05 second steps
  37. dt = 0.05
  38. t = np.arange(0.0, 20, dt)
  40. # th1 and th2 are the initial angles (degrees)
  41. # w10 and w20 are the initial angular velocities (degrees per second)
  42. th1 = 120.0
  43. w1 = 0.0
  44. th2 = -10.0
  45. w2 = 0.0
  47. # initial state
  48. state = np.radians([th1, w1, th2, w2])
  50. # integrate your ODE using scipy.integrate.
  51. y = integrate.odeint(derivs, state, t)
  53. x1 = L1*sin(y[:, 0])
  54. y1 = -L1*cos(y[:, 0])
  56. x2 = L2*sin(y[:, 2]) + x1
  57. y2 = -L2*cos(y[:, 2]) + y1
  59. fig = plt.figure()
  60. ax = fig.add_subplot(111, autoscale_on=False, xlim=(-2, 2), ylim=(-2, 2))
  61. ax.set_aspect('equal')
  62. ax.grid()
  64. line, = ax.plot([], [], 'o-', lw=2)
  65. time_template = 'time = %.1fs'
  66. time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
  69. def init():
  70.     line.set_data([], [])
  71.     time_text.set_text('')
  72.     return line, time_text
  75. def animate(i):
  76.     thisx = [0, x1[i], x2[i]]
  77.     thisy = [0, y1[i], y2[i]]
  79.     line.set_data(thisx, thisy)
  80.     time_text.set_text(time_template % (i*dt))
  81.     return line, time_text
  83. ani = animation.FuncAnimation(fig, animate, np.arange(1, len(y)),
  84.                               interval=25, blit=True, init_func=init)
  86. plt.show()

下载这个示例