Trouble animating a double pendulum

Hello! For a project in my physics class, I am animating a double pendulum, but I keep having trouble with the animation. It keeps accelerating and it doesn't seem accurate. I showed the code to my physics teacher, and he said that it all seemed fine and couldn't understand why it isn't working. Should I stick with the method I am using right now or switch to using ode45?

syms a_1 a_2

theta_1 = input('What do you want Theta 1 to equal? ')*pi/180;
theta_2 = input('What do you want Theta 2 to equal? ')*pi/180;

l_1 = 1;
l_2 = 1;

g = 9.8;

m_1 = 1;
m_2 = 1;

dt = .01;

w_1 = 0;
w_2 = 0;

for n=1:1000
   clf
   a_1=-(g*(2*m_1+m_2)*sin(theta_1)+m_2*g*sin(theta_1-2*theta_2)+2*sin(theta_1-theta_2)*m_2*(((w_2)^2)*l_2+((w_1)^2)*l_1*cos(theta_1-theta_2)))/(l_1*2*m_1+m_2-m_2*cos(2*(theta_1-theta_2)));
   a_2 = (2*sin(theta_1-theta_2)*((w_1)^2)*l_1*(m_1+m_2)+g*(m_1+m_2)*cos(theta_1)+((w_2)^2)*l_2*m_2*cos(theta_1-theta_2))/(l_2*(2*m_1+m_2-m_2*cos(2*theta_1-2*theta_2)));
   
   w_1 = w_1+a_1*dt;
   w_2 = w_2+a_2*dt;

   theta_1 = theta_1+w_1*dt;
   theta_2 = theta_2+w_2*dt;

   x_1 = l_1*sin(theta_1);
   x_2 = l_1*sin(theta_1)+l_2*sin(theta_2);

   y_1 = -l_1*cos(theta_1);
   y_2 = -l_1*cos(theta_1)-l_2*cos(theta_2);

   plot(0, 0, 'O')
   hold on
   plot(x_1, y_1, '*')
   hold on
   line([0 x_1], [0 y_1])
   hold on
   plot(x_2, y_2, 'o')
   hold on
   line([x_1 x_2], [y_1 y_2])
   axis([-10 10 -10 10])
   pause(.01)
end

[Torsten]

Should I stick with the method I am using right now or switch to using ode45?

If you solve a problem with a self-written numerical code, you should always compare with a professional solver if possible.

Concerning your code:

Don't define a1 and a2 as symbolic.

And don't plot within the loop. Save x_1, x_2, y_1 and y_2 in an array of length 1001 and plot after the loop.

Answer 1

[Torsten]

theta_1 = 19*pi/180;
theta_2 = 30*pi/180;

l_1 = 1;
l_2 = 1;

g = 9.8;

m_1 = 1;
m_2 = 1;

M = m_1 + m_2;

dt = .001;

w_1 = 0;
w_2 = 0;

x_1 = zeros(10000,1);
y_1 = zeros(10000,1);
x_2 = zeros(10000,1);
y_2 = zeros(10000,1);

x_1(1) = l_1*sin(theta_1);
x_2(1) = x_1(1) + l_2*sin(theta_2);

y_1(1) = -l_1*cos(theta_1);
y_2(1) = y_1(1) - l_2*cos(theta_2);

for n=2:10000

   % Source: https://physics.umd.edu/hep/drew/pendulum2.html
   dtheta = theta_1 - theta_2;
   alpha = m_1 + m_2*sin(dtheta)^2;
   a_1 = (-sin(dtheta)*(m_2*l_1*w_1^2*cos(dtheta)+m_2*l_2*w_2^2)-...
          g*(M*sin(theta_1)-m_2*sin(theta_2)*cos(dtheta)))/(l_1*alpha);
   a_2 = (sin(dtheta)*(M*l_1*w_1^2+m_2*l_2*w_2^2*cos(dtheta))+...
          g*M*(sin(theta_1)*cos(dtheta)-sin(theta_2)))/(l_2*alpha);
   %dtheta = theta_1 - theta_2;
   %Mat = [M*l_1,m_2*l_2*cos(dtheta);l_1*cos(dtheta),l_2];
   %rhs = [-m_2*l_2*w_2^2*sin(dtheta)-M*g*sin(theta_1);l_1*w_1^2*sin(dtheta)-g*sin(theta_2)];
   %sol = Mat\rhs;
   %a_1 = sol(1);
   %a_2 = sol(2);
           
   w_1_new = w_1+a_1*dt;
   w_2_new = w_2+a_2*dt;

   theta_1_new = theta_1+w_1*dt;
   theta_2_new = theta_2+w_2*dt;

   w_1 = w_1_new;
   w_2 = w_2_new;
   
   theta_1 = theta_1_new;
   theta_2 = theta_2_new;
   
   x_1(n) = l_1*sin(theta_1);
   x_2(n) = x_1(n) + l_2*sin(theta_2);

   y_1(n) = -l_1*cos(theta_1);
   y_2(n) = y_1(n) - l_2*cos(theta_2);


end
hold on
plot(x_1,y_1)
plot(x_2,y_2)
hold off
grid on

[Sam Chak]

A double pendulum is a chaotic system, and the second pendulum will exhibit unpredictable patterns of motion. If the simulation responses appear predictable, either the equations are incorrect or the ratios of lengths and masses are carefully tuned, resulting in more predictable behavior.

Chaotic behavior:

Predictable behavior: