77import scipy .integrate as integrate
88import matplotlib .animation as animation
99
10- G = 9.8 # acceleration due to gravity, in m/s^2
11- L1 = 1.0 # length of pendulum 1 in m
12- L2 = 1.0 # length of pendulum 2 in m
13- M1 = 1.0 # mass of pendulum 1 in kg
14- M2 = 1.0 # mass of pendulum 2 in kg
10+ G = 9.8 # acceleration due to gravity, in m/s^2
11+ L1 = 1.0 # length of pendulum 1 in m
12+ L2 = 1.0 # length of pendulum 2 in m
13+ M1 = 1.0 # mass of pendulum 1 in kg
14+ M2 = 1.0 # mass of pendulum 2 in kg
1515
1616
1717def derivs (state , t ):
1818
1919 dydx = np .zeros_like (state )
2020 dydx [0 ] = state [1 ]
2121
22- del_ = state [2 ] - state [0 ]
23- den1 = (M1 + M2 ) * L1 - M2 * L1 * cos (del_ ) * cos (del_ )
24- dydx [1 ] = (M2 * L1 * state [1 ] * state [1 ] * sin (del_ ) * cos (del_ )
25- + M2 * G * sin (state [2 ]) * cos (del_ ) + M2 *
26- L2 * state [3 ] * state [3 ] * sin (del_ )
27- - (M1 + M2 ) * G * sin (state [0 ])) / den1
22+ del_ = state [2 ]- state [0 ]
23+ den1 = (M1 + M2 )* L1 - M2 * L1 * cos (del_ )* cos (del_ )
24+ dydx [1 ] = (M2 * L1 * state [1 ]* state [1 ]* sin (del_ )* cos (del_ )
25+ + M2 * G * sin (state [2 ])* cos (del_ ) + M2 * L2 * state [3 ]* state [3 ]* sin (del_ )
26+ - (M1 + M2 )* G * sin (state [0 ]))/ den1
2827
2928 dydx [2 ] = state [3 ]
3029
31- den2 = (L2 / L1 ) * den1
32- dydx [3 ] = (- M2 * L2 * state [3 ] * state [3 ] * sin (del_ ) * cos (del_ )
33- + (M1 + M2 ) * G * sin (state [0 ]) * cos (del_ )
34- - (M1 + M2 ) * L1 * state [1 ] * state [1 ] * sin (del_ )
35- - (M1 + M2 ) * G * sin (state [2 ])) / den2
30+ den2 = (L2 / L1 )* den1
31+ dydx [3 ] = (- M2 * L2 * state [3 ]* state [3 ]* sin (del_ )* cos (del_ )
32+ + (M1 + M2 )* G * sin (state [0 ])* cos (del_ )
33+ - (M1 + M2 )* L1 * state [1 ]* state [1 ]* sin (del_ )
34+ - (M1 + M2 )* G * sin (state [2 ]))/ den2
3635
3736 return dydx
3837
@@ -47,17 +46,19 @@ def derivs(state, t):
4746th2 = - 10.0
4847w2 = 0.0
4948
49+ rad = pi / 180
50+
5051# initial state
51- state = np .radians ([th1 , w1 , th2 , w2 ])
52+ state = np .array ([th1 , w1 , th2 , w2 ])* pi / 180.
5253
5354# integrate your ODE using scipy.integrate.
5455y = integrate .odeint (derivs , state , t )
5556
56- x1 = L1 * sin (y [:, 0 ])
57- y1 = - L1 * cos (y [:, 0 ])
57+ x1 = L1 * sin (y [:,0 ])
58+ y1 = - L1 * cos (y [:,0 ])
5859
59- x2 = L2 * sin (y [:, 2 ]) + x1
60- y2 = - L2 * cos (y [:, 2 ]) + y1
60+ x2 = L2 * sin (y [:,2 ]) + x1
61+ y2 = - L2 * cos (y [:,2 ]) + y1
6162
6263fig = plt .figure ()
6364ax = fig .add_subplot (111 , autoscale_on = False , xlim = (- 2 , 2 ), ylim = (- 2 , 2 ))
@@ -67,23 +68,21 @@ def derivs(state, t):
6768time_template = 'time = %.1fs'
6869time_text = ax .text (0.05 , 0.9 , '' , transform = ax .transAxes )
6970
70-
7171def init ():
7272 line .set_data ([], [])
7373 time_text .set_text ('' )
7474 return line , time_text
7575
76-
7776def animate (i ):
7877 thisx = [0 , x1 [i ], x2 [i ]]
7978 thisy = [0 , y1 [i ], y2 [i ]]
8079
8180 line .set_data (thisx , thisy )
82- time_text .set_text (time_template % ( i * dt ))
81+ time_text .set_text (time_template % ( i * dt ))
8382 return line , time_text
8483
8584ani = animation .FuncAnimation (fig , animate , np .arange (1 , len (y )),
86- interval = 25 , blit = True , init_func = init )
85+ interval = 25 , blit = True , init_func = init )
8786
8887#ani.save('double_pendulum.mp4', fps=15)
8988plt .show ()
0 commit comments