0

I am using the following code to solve the system of differential equations numerically and then animate it:

for i in range(steps-1):
    Theta1 = t1[i]
    Theta2 = t2[i]
    dTheta1 = w1[i]
    dTheta2 = w2[i]
    a1 = (g*(np.sin(Theta2)*np.cos(Theta1-Theta2)-mu*np.sin(Theta1))-(l2*dTheta2*dTheta2+l1*dTheta1*dTheta1*np.cos(Theta1-Theta2))*np.sin(Theta1-Theta2))/(l1*(mu-np.cos(Theta1-Theta2)*np.cos(Theta1-Theta2)))
    a2 = (mu*g*(np.sin(Theta1)*np.cos(Theta1-Theta2)-np.sin(Theta2))+(mu*l1*dTheta1*dTheta1+l2*dTheta2*dTheta2*np.cos(Theta1-Theta2))*np.sin(Theta1-Theta2))/(l2*(mu-np.cos(Theta1-Theta2)*np.cos(Theta1-Theta2)))
    w1[i+1] = w1[i] + dt*a1
    w2[i+1] = w2[i] + dt*a2
    t1[i+1] = t1[i] + dt*w1[i]

    t2[i+1] = t2[i] + dt*w2[i]

This gives me errors such as RuntimeWarning: invalid value encountered in longdouble_scalars and RuntimeWarning: overflow encountered in longdouble_scalars. At first, I thought this might be happening due to precisions errors, so I tried using longdouble. But that didn't help at all. When I run the animation, I get a few frames of unphysical systems and then it shuts downs. I have checked my equation from multiple sources and I ended up using their equations too and the problem still persists. What am I doing wrong?

Edit: Okay so simply decreasing the step size seemed to not give me the error any more, but now animation is unphysical. The pendulum seems to pick up speed over time, and also performs full rotations, even when it started from rest at an extreme. I am using a 4th order RK method and still getting this problem. Any ideas? I have provided my new code below

def diff_eq_a1(theta1,theta2,w1,w2):
    a11 = mu
    a12 = np.cos(theta1-theta2)
    a21 = np.cos(theta1-theta2)
    a22 = 1
    
    b1 = (mu*g*l1*np.sin(theta1)) + (l2*np.sin(theta1-theta2)*w2*w2)
    b2 = (g*l2*np.sin(theta2))-(l1*np.sin(theta1-theta2)*w1*w1)
    
    return -(1/l1)*(b1*a22-b2*a12)/(a11*a22-a21*a12)
    

def diff_eq_a2(theta1,theta2,w1,w2):
    a11 = mu
    a12 = np.cos(theta1-theta2)
    a21 = np.cos(theta1-theta2)
    a22 = 1
    
    b1 = (mu*g*l1*np.sin(theta1)) + (l2*np.sin(theta1-theta2)*w2*w2)
    b2 = (g*l2*np.sin(theta2))-(l1*np.sin(theta1-theta2)*w1*w1)
    
    return -(1/l2)*(b2*a11-b1*a21)/(a11*a22-a21*a12)


for i in range(steps-1):
    a1_k1 = dt*diff_eq_a1(t1[i],t2[i],w1[i],w2[i])
    a1_k2 = dt*diff_eq_a1((t1[i]+0.5*a1_k1),t2[i],w1[i],w2[i])
    a1_k3 = dt*diff_eq_a1((t1[i]+0.5*a1_k2),t2[i],w1[i],w2[i])
    a1_k4 = dt*diff_eq_a1((t1[i]+a1_k3),t2[i],w1[i],w2[i])
    
    w1[i+1] = w1[i] + (a1_k1/6)+(a1_k2/3)+(a1_k3/3)+(a1_k4/6)
    
    a2_k1 = dt*diff_eq_a1(t1[i],t2[i],w1[i],w2[i])
    a2_k2 = dt*diff_eq_a1(t1[i],(t2[i]+0.5*a2_k1),w1[i],w2[i])
    a2_k3 = dt*diff_eq_a1(t1[i],(t2[i]+0.5*a2_k2),w1[i],w2[i])
    a2_k4 = dt*diff_eq_a1(t1[i],(t2[i]+a2_k3),w1[i],w2[i])
    
    w2[i+1] = w2[i] + (a2_k1/6)+(a2_k2/3)+(a2_k3/3)+(a2_k4/6)
    
    t1[i+1] = t1[i] + dt*w1[i]
    t2[i+1] = t2[i] + dt*w2[i]

x1 = l1*np.sin(t1)
x2 = x1 + l2*np.sin(t2)
y1 = - l1*np.cos(t1)
y2 = y1 - l2*np.cos(t2)

I know its not very neat right now, but I really just want to get it working first

Souroy
  • 21
  • 2
  • What is your step size? The Euler method tends to leave on a tangent, which for the sharp turns of a double pendulum can mean a rapid growth in errors, also in unphysical directions. Try a higher order method. Use a method with local-error estimates to select an appropriate step size for the sharp turns. – Lutz Lehmann Jan 13 '21 at 17:49
  • See https://stackoverflow.com/questions/65224923/i-want-to-have-the-pendulum-blob-in-my-double-pendulum for a successful implementation, also https://math.stackexchange.com/questions/2084922/transform-the-double-pendulum-differential-equations-into-a-first-order for a better organization of the long terms – Lutz Lehmann Jan 13 '21 at 17:56
  • Combine both acceleration computations in one function, then apply RK4 for a mechanical system like in https://stackoverflow.com/questions/60405185/is-there-a-better-way-to-tidy-this-chuck-of-code-it-is-the-runge-kutta-4th-orde, or in without using vectors https://stackoverflow.com/questions/57253590/plot-orbit-of-two-body-using-rk4/ – Lutz Lehmann Jan 14 '21 at 08:14
  • Why your approach does not work (... much better than Euler): https://stackoverflow.com/a/62625121/3088138, https://stackoverflow.com/a/53996811/3088138, https://stackoverflow.com/a/58258912/3088138 – Lutz Lehmann Jan 14 '21 at 08:37

0 Answers0