7
7
import scipy .integrate as integrate
8
8
import matplotlib .animation as animation
9
9
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
15
15
16
16
17
17
def derivs (state , t ):
18
18
19
19
dydx = np .zeros_like (state )
20
20
dydx [0 ] = state [1 ]
21
21
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
27
28
28
29
dydx [2 ] = state [3 ]
29
30
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
35
36
36
37
return dydx
37
38
@@ -46,19 +47,17 @@ def derivs(state, t):
46
47
th2 = - 10.0
47
48
w2 = 0.0
48
49
49
- rad = pi / 180
50
-
51
50
# initial state
52
- state = np .array ([th1 , w1 , th2 , w2 ])* pi / 180.
51
+ state = np .radians ([th1 , w1 , th2 , w2 ])
53
52
54
53
# integrate your ODE using scipy.integrate.
55
54
y = integrate .odeint (derivs , state , t )
56
55
57
- x1 = L1 * sin (y [:,0 ])
58
- y1 = - L1 * cos (y [:,0 ])
56
+ x1 = L1 * sin (y [:, 0 ])
57
+ y1 = - L1 * cos (y [:, 0 ])
59
58
60
- x2 = L2 * sin (y [:,2 ]) + x1
61
- y2 = - L2 * cos (y [:,2 ]) + y1
59
+ x2 = L2 * sin (y [:, 2 ]) + x1
60
+ y2 = - L2 * cos (y [:, 2 ]) + y1
62
61
63
62
fig = plt .figure ()
64
63
ax = fig .add_subplot (111 , autoscale_on = False , xlim = (- 2 , 2 ), ylim = (- 2 , 2 ))
@@ -68,21 +67,23 @@ def derivs(state, t):
68
67
time_template = 'time = %.1fs'
69
68
time_text = ax .text (0.05 , 0.9 , '' , transform = ax .transAxes )
70
69
70
+
71
71
def init ():
72
72
line .set_data ([], [])
73
73
time_text .set_text ('' )
74
74
return line , time_text
75
75
76
+
76
77
def animate (i ):
77
78
thisx = [0 , x1 [i ], x2 [i ]]
78
79
thisy = [0 , y1 [i ], y2 [i ]]
79
80
80
81
line .set_data (thisx , thisy )
81
- time_text .set_text (time_template % ( i * dt ))
82
+ time_text .set_text (time_template % ( i * dt ))
82
83
return line , time_text
83
84
84
85
ani = animation .FuncAnimation (fig , animate , np .arange (1 , len (y )),
85
- interval = 25 , blit = True , init_func = init )
86
+ interval = 25 , blit = True , init_func = init )
86
87
87
88
#ani.save('double_pendulum.mp4', fps=15)
88
89
plt .show ()
0 commit comments