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 * L2 * state [3 ]* state [3 ]* sin (del_ )
26- - (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 *
26+ L2 * state [3 ] * state [3 ] * sin (del_ )
27+ - (M1 + M2 ) * G * sin (state [0 ])) / den1
2728
2829 dydx [2 ] = state [3 ]
2930
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
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
3536
3637 return dydx
3738
@@ -46,19 +47,19 @@ def derivs(state, t):
4647th2 = - 10.0
4748w2 = 0.0
4849
49- rad = pi / 180
50+ rad = pi / 180
5051
5152# initial state
52- state = np .array ([th1 , w1 , th2 , w2 ])* pi / 180.
53+ state = np .array ([th1 , w1 , th2 , w2 ]) * pi / 180.
5354
5455# integrate your ODE using scipy.integrate.
5556y = integrate .odeint (derivs , state , t )
5657
57- x1 = L1 * sin (y [:,0 ])
58- y1 = - L1 * cos (y [:,0 ])
58+ x1 = L1 * sin (y [:, 0 ])
59+ y1 = - L1 * cos (y [:, 0 ])
5960
60- x2 = L2 * sin (y [:,2 ]) + x1
61- y2 = - L2 * cos (y [:,2 ]) + y1
61+ x2 = L2 * sin (y [:, 2 ]) + x1
62+ y2 = - L2 * cos (y [:, 2 ]) + y1
6263
6364fig = plt .figure ()
6465ax = fig .add_subplot (111 , autoscale_on = False , xlim = (- 2 , 2 ), ylim = (- 2 , 2 ))
@@ -68,21 +69,23 @@ def derivs(state, t):
6869time_template = 'time = %.1fs'
6970time_text = ax .text (0.05 , 0.9 , '' , transform = ax .transAxes )
7071
72+
7173def init ():
7274 line .set_data ([], [])
7375 time_text .set_text ('' )
7476 return line , time_text
7577
78+
7679def animate (i ):
7780 thisx = [0 , x1 [i ], x2 [i ]]
7881 thisy = [0 , y1 [i ], y2 [i ]]
7982
8083 line .set_data (thisx , thisy )
81- time_text .set_text (time_template % ( i * dt ))
84+ time_text .set_text (time_template % ( i * dt ))
8285 return line , time_text
8386
8487ani = animation .FuncAnimation (fig , animate , np .arange (1 , len (y )),
85- interval = 25 , blit = True , init_func = init )
88+ interval = 25 , blit = True , init_func = init )
8689
8790#ani.save('double_pendulum.mp4', fps=15)
8891plt .show ()
0 commit comments