In [1]:

    

using PyPlot

g = 9.7864
m1 = 1
m2 = 0.05
l1 = 0.5
l2 = 0.2


    

    
Out[1]:





0.2

In [2]:

    

der_ang(p1, p2, c1, c2) = p1 * c1 + p2 * c2
der_p(ang, c1, c2) = c1 * sin(ang) + c2
fixang(ang) = atan2(sin(ang),cos(ang))


    

    
Out[2]:





fixang (generic function with 1 method)

In [24]:

    

dt = 0.0001;
t_max = 500;
count = Int64(t_max/dt);
x=0

#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 = 1 : 8
#condicoes iniciais envolvendo energia
E_0 = 2 + (j-1) * 4;
ang1 = 0;
    ang1l = 0
ang2 = pi/2;
    ang2l = 0
p1 = sqrt(2 * m1 * l1^2 *(E_0 - m2 * l2 * g));
    p1l = 0
p2 = 0;  
    p2l = 0
print(E_0);
    print("\n")
for i = 1 : count-1
cv1 = -cos(ang1 - ang2) * cp;
cv2 = -sin(ang1 - ang2) * cp * p1 * p2;
k_ang1_1 = der_ang(p1, p2, c1, cv1) * dt;
k_ang2_1 = der_ang(p2, p1, c2, cv1) * dt;
k_p1_1 = der_p(ang1, c3, cv2) * dt;
k_p2_1 = der_p(ang2, c4, -cv2) * dt;
xang1 = ang1 + k_ang1_1;
xang2 = ang2 + k_ang2_1;
xp1 = p1 + k_p1_1;
xp2 = p2 + 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;
        ang1l = ang1
        ang2l = ang2
ang1 = fixang(ang1 + k_ang1_2);
ang2 = fixang(ang2 + k_ang2_2);
        p1 = p1 + k_p1_2;
p2 = p2 + k_p2_2;
        if( sign(ang2) != sign(ang2l))
if( p2 > 0)
x = x + 1;
                push!(out_ang1, ang1)
                push!(out_p1, p1)
end
end

end
end


    

    




2
4
6
8
10
12
14
16
18

In [26]:

    

plot(out_ang1/pi, out_p1, ".k", markersize=1)
xlabel("p_1")
ylabel("Omega_1 / pi")


    

    













    
Out[26]:





PyObject <matplotlib.text.Text object at 0x7f80dde8fa50>

In [ ]: