-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path2_double_pendulums.py
130 lines (105 loc) · 5.44 KB
/
2_double_pendulums.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from matplotlib.animation import FuncAnimation
from matplotlib import rcParams
def derivatives(y, t, L1, L2, m1, m2):
theta1, omega1, theta2, omega2 = y
g = 9.81
dtheta1_dt = omega1
domega1_dt = (-g*(2*m1 + m2)*np.sin(theta1) - m2*g*np.sin(theta1 - 2*theta2) - 2*np.sin(theta1 - theta2)*m2*(omega2**2*L2 + omega1**2*L1*np.cos(theta1 - theta2)))/(L1*(2*m1 + m2 - m2*np.cos(2*theta1 - 2*theta2)))
dtheta2_dt = omega2
domega2_dt = (2*np.sin(theta1 - theta2)*(omega1**2*L1*(m1 + m2) + g*(m1 + m2)*np.cos(theta1) + omega2**2*L2*m2*np.cos(theta1 - theta2)))/(L2*(2*m1 + m2 - m2*np.cos(2*theta1 - 2*theta2)))
return dtheta1_dt, domega1_dt, dtheta2_dt, domega2_dt
#initial conditions
theta1 = np.radians(90)
theta2 = np.radians(45)
omega1 = 0
omega2 = 0
theta3 = np.radians(100)
theta4 = np.radians(55)
omega3 = 0
omega4 = 0
y0 = [theta1, theta2, omega1, omega2]
y1 = [theta3, theta4, omega3, omega4]
#odient
t = np.linspace(0,10,30*10)
# parameters for the pendulums
L1 = 0.8 # length of pendulum 1 in meters
L2 = 0.8 # length of pendulum 2 in meters
m1 = 1 # mass of pendulum 1 in kilograms
m2 = 1 # mass of pendulum 2 in kilograms
#parameters for 2nd pendulum
L3 = 0.8 # length of pendulum 1 in meters
L4 = 0.8 # length of pendulum 2 in meters
m3 = 1 # mass of pendulum 1 in kilograms
m4 = 1 # mass of pendulum 2 in kilograms
# solve the equations of motion
solution = odeint(derivatives, y0, t, args=(L1, L2, m1, m2))
solution2 = odeint(derivatives, y1, t, args=(L3, L4, m3, m4))
theta1 = solution[:,0]
theta2 = solution[:,2]
theta3 = solution2[:,0]
theta4 = solution2[:,2]
#x,y coordinates
x1 = L1*np.sin(theta1)
y1 = -L1*np.cos(theta1)
x2 = L1*np.sin(theta1) + L2*np.sin(theta2)
y2 = -(L1*np.cos(theta1) + L2*np.cos(theta2))
x3 = L3*np.sin(theta3)
y3 = -L3*np.cos(theta3)
x4 = L3*np.sin(theta3) + L4*np.sin(theta4)
y4 = -(L3*np.cos(theta3) + L4*np.cos(theta4))
#subplot
def init():
fig, ax = plt.subplots() # create a figure and a subplot
ax.set_xlim(-L1 - L2, L1 + L2) # set the limits of the x axis
ax.set_ylim(-0.2, L1 + L2) # set the limits of the y axis
ax.set_xticks(np.arange(-2, 2, 0.3)) # set the x ticks
ax.set_yticks(np.arange(-2, 2, 0.3)) # set the y ticks
ax.grid(True, color='gray', linewidth=0.5) # add a grid with gray color and thinner lines
ax.set_facecolor('black') # set the background color to black
ax.tick_params(axis='both', colors='white') # set the color of the tick labels to white
ax.plot(0, 0, 'ko', markersize = 10, color = 'silver') # add a black circle at the origin
# plot the arms of the pendulums with thicker lines and different colors
line1, = ax.plot([0, x1[0]], [0, y1[0]], 'darkblue', linewidth = 2) # 'darkblue' for the pendulum arm
line2, = ax.plot([x1[0], x2[0]], [y1[0], y2[0]], 'darkorange', linewidth = 2) # 'darkorange' for the pendulum arm
line3, = ax.plot([0, x3[0]], [0, y3[0]], 'blue', linewidth = 2) # 'darkblue' for the pendulum arm
line4, = ax.plot([x3[0], x4[0]], [y3[0], y4[0]], 'orange', linewidth = 2) # 'darkorange' for the pendulum arm
# plot the mass blobs
blob1, = ax.plot([x1[0]], [y1[0]], 'o', markersize = m1*8 ,color='darkblue') # 'darkblue' for the blob
blob2, = ax.plot([x2[0]], [y2[0]], 'o', markersize = m2*10,color='darkorange') # 'darkorange' for the blob
blob3, = ax.plot([x3[0]], [y3[0]], 'o', markersize = m3*8 ,color='blue') # 'darkblue' for the blob
blob4, = ax.plot([x4[0]], [y4[0]], 'o', markersize = m4*10,color='orange') # 'darkorange' for the blob
#trails
#trail1, = ax.plot([], [], 'c-', linewidth=1) # cyan trail for bob 1
trail2, = ax.plot([], [], 'm-', linewidth=1) # magenta trail for bob 2
#trail3, = ax.plot([], [], 'b-', linewidth=1) # cyan trail for bob 1
trail4, = ax.plot([], [], 'y-', linewidth=1) # magenta trail for bob 2
return fig, line1, line2, line3, line4, blob1, blob2, blob3, blob4, trail2, trail4
def update(frame):
line1.set_data([0, x1[frame]], [0, y1[frame]])
line2.set_data([x1[frame], x2[frame]], [y1[frame], y2[frame]])
line3.set_data([0, x3[frame]], [0, y3[frame]])
line4.set_data([x3[frame], x4[frame]], [y3[frame], y4[frame]])
blob1.set_data(x1[frame], y1[frame])
blob2.set_data(x2[frame], y2[frame])
blob3.set_data(x3[frame], y3[frame])
blob4.set_data(x4[frame], y4[frame])
# update the data for the trails
# limit the length of the trails
trail_length = 10 # adjust this to change the length of the trails
if frame > trail_length:
#trail1.set_data(x1[frame-trail_length:frame], y1[frame-trail_length:frame])
trail2.set_data(x2[frame-trail_length:frame], y2[frame-trail_length:frame])
#trail3.set_data(x3[frame-trail_length:frame], y3[frame-trail_length:frame])
trail4.set_data(x4[frame-trail_length:frame], y4[frame-trail_length:frame])
else:
#trail1.set_data(x1[:frame], y1[:frame])
trail2.set_data(x2[:frame], y2[:frame])
#trail3.set_data(x3[:frame], y3[:frame])
trail4.set_data(x4[:frame], y4[:frame])
return fig, line1, line2, line3, line4, blob1, blob2, blob3, blob4, trail2, trail4
fig, line1, line2, line3, line4, blob1, blob2, blob3, blob4, trail2, trail4 = init()
animation = FuncAnimation(fig, update, frames=np.arange(len(t)), init_func=init, interval=200, blit=True)
animation.save('ml_work/double_pendulum/outputs/2double_pendulum13.gif', writer='imagemagick')