In [50]:
import numpy as np
import matplotlib.pyplot as plt
g = 9.7864
m1 = 1
m2 = 0.05
l1 = 0.5
l2 = 0.2
In [51]:
def der_ang(p1, p2, c1, c2):
return p1 * c1 + p2 * c2
def der_p(ang, c1, c2):
return c1 * np.sin(ang) + c2
def fixang(ang):
return np.arctan2(np.sin(ang),np.cos(ang))
In [96]:
dt = 0.01;
t_max = 100;
count = t_max/dt;
x=0
ang1 = np.zeros((count, 9));
ang2 = np.zeros((count, 9));
p1 = np.zeros((count, 9));
p2 = np.zeros((count, 9));
#constantes de uso comum
c1 = 1 / (m1 * l1 ** 2);
c2 = 1 / (m2 * l2 ** 2);
c3 = -1 * m1 * g * l1;
c4 = -1 * m2 * g * l2;
cp = 1 / (m1 * l1 * l2);
out_p1 = [];
out_ang1 =[];
for j in range(1, 9):
#condicoes iniciais envolvendo energia
E_0 = 0.2 + (j-1) * 0.2;
ang1[1, j] = 0;
ang2[1, j] = np.pi/2;
p1[1, j] = np.sqrt(2 * m1 * l1 ** 2 *(E_0 - m2 * l2 * g));
p2[1, j] = 0;
print(E_0);
for i in range(1, int(count-1)):
cv1 = -np.cos(ang1[i, j] - ang2[i, j]) * cp;
cv2 = -np.sin(ang1[i, j] - ang2[i, j]) * cp * p1[i, j] * p2[i, j];
k_ang1_1 = der_ang(p1[i, j], p2[i, j], c1, cv1) * dt;
k_ang2_1 = der_ang(p2[i, j], p1[i, j], c2, cv1) * dt;
k_p1_1 = der_p(ang1[i, j], c3, cv2) * dt;
k_p2_1 = der_p(ang2[i, j], c4, -cv2) * dt;
xang1 = ang1[i, j] + k_ang1_1;
xang2 = ang2[i, j] + k_ang2_1;
xp1 = p1[i, j] + k_p1_1;
xp2 = p2[i, j] + k_p2_1;
k_ang1_2 = der_ang(xp1, xp2, c1, cv1) * dt;
k_ang2_2 = der_ang(xp2, xp1, c2, cv1) * dt;
k_p1_2 = der_p(xang1, c3, cv2) * dt;
k_p2_2 = der_p(xang2, c4, -cv2) * dt;
ang1[i+1, j] = fixang(ang1[i, j] + k_ang1_2);
ang2[i+1, j] = fixang(ang2[i, j] + k_ang2_2);
if( np.sign(ang1[i+1, j]) != np.sign(ang1[i, j])):
if( p2[i, j] > 0):
x = x + 1;
out_ang1.append(ang2[i, j])
out_p1.append(p2[i, j])
p1[i+1, j] = p1[i, j] + k_p1_2;
p2[i+1, j] = p2[i, j] + k_p2_2;
In [97]:
plt.plot(np.array(out_ang1) / np.pi, out_p1, '.k')
plt.show()
In [92]:
plt.plot(ang2, p2, '.k')
plt.show()
In [ ]: