How to run the animation of double pendulum chaotic nature?

I was writing a code for simulating the chaotic nature of double pendulum using runge-kutta 4th order method. First of all,i wrote the code for the calculation of the k's. It worked but then i wrote the code for the animation part but the animation is not working right. Now,i am starting to doubt if the calculation part is right. Anyone,who could check my code and tell me where did i go wrong?

function double_pendulum(m1,m2,l1,l2,theta1,theta2,tmax)
%value of the constants
  m1=0.1;
  m2=0.1;
  l1=0.2;
  l2=0.3;
  theta1=30;
  theta2=40;
  g=9.8;
  tmax=10;
 %RK matrices
  c=[0 1/2 1/2 1];
  a=[0 0 0 0;0 1/2 0 0;0 0 1/2 0;0 0 0 1];
  w=[1/6;2/6;2/6;1/6];
  theta1=theta1*pi/180;
  theta2=theta2*pi/180;
  s(:,1)=[theta1 theta2 0 0]';
  F=@(t,s)[s(3) s(4) -m2*l1*(s(4))^2*sin(s(1)-s(2))*cos(s(1)-s(2))+g*m2*sin(s(2))*cos(s(1)-s(2))-m2*l2*(s(4)^2)*sin(s(1)-s(2))-(m1+m2)*g*sin(s(1))/l1*(m1+m2)-m2*l1*(cos(s(1)-s(2))^2) m2*l2*(s(4)^2)*sin(s(1)-s(2))*cos(s(1)-s(2))+g*sin(s(1))*cos(s(1)-s(2))*(m1+m2)+l1*(s(3))^2*sin(s(1)-s(2))*(m1+m2)-g*sin(s(2))*(m1+m2)/l2*(m1+m2)-m2*l2*cos(s(1)-s(2))^2]';
  h=0.1;
  i=1;
  t(1)=0;
%==================================================
%Computation of position and angular Velocity
while t(i)<tmax 
  k=zeros (length(s(:,1)),length(c));
  for j=1 : length(c)

        k(:,j)=h*F(t(i)+h*c(j),s(:,i)+k*a(j,:)')

    end
     s(:,i+1)=s(:,i)+k*w;
     t(i+1)=t(i)+h; 
     i=i+1;
end
    x1=l1*sin(s(1,:));
    x2=l1*sin(s(1,:))+l2*sin(s(2,:));      % x- axis position of pendulum...... s(1,:)= All theta values
    y1=-l1*cos(s(1,:));  %y-axis position of pendulum
    y2=-l1*cos(s(1,:))-l2*cos(s(2,:));
      for i=1:length(s(1,:))
        for j=1:length(s(2,:))
        clf;
        plot(x1(i),y1(i),'o','markersize',30,'markerfacecolor','m');  % plots the bob1
        plot([0 x1(i)],[l1 y1(i)],'linewidth',3,'color','b');    %plots rod1
        plot(x2(j),y2(j),'o','markersize',30,'markerfacecolor','m');%plots bob2
        plot([x1(i) x2(j)],[y1(i) y2(j) ],'linewidth',3,'color','r'); %plots rod2
        axis([-2*l1 2*l1 -1 2*l2]);
        hold on;                                                                                %and other value as per x(i)...so while
                                                                                             %ploting first point is (0,L) and joins
        axis off;                                                                             %this point with other point.
        plot([-5 5],[l1 l1],'k','linewidth',3);                                               %plots the top line
        hold off;
        pause(0.1);
        end 
      end

Answer 1

[Jan]

Compare the results with the ODE45 integrator:

[T, Y] = ode45(F, [0, tmax], s(:, 1).');
figure;
plot(T, Y);
figure;
plot(s.')  % Your solution

The trajectories are similar, but not equal. This is a problem of the too large step width h=0.1. As soon as you reduce the step width to h=0.01 the results become nearly equal. Your Runge Kutta seems to be fine.

while t < tmax create the last time point 10.1, not 10.0. For a constant step width prefer:

t = linspace(0, 3, round(3 / 0.1) + 1);
s = zeros(4, numel(t));   % Pre-allocate!
for i = 2:numel(t)
  ...
end

Letting an array grow iteratively requires a lot of resources nd can reduce the runtime dramatically. Therefore a pre-allocation of the output is a good programming practice.

Start the loop at time point 2, because at t(1) the initial values are defined already. The wrong last time point is a small inaccuracy, but it can be avoided easily.

But you are right: I do not trust the trajectories also. If the integration is okay, the ODE must have a problem. How did you obtain it? Compare it with e.g. http://sophia.dtp.fmph.uniba.sk/~kovacik/doublePendulum.pdf. E.g. in the 4th component you have:

 ... - g*sin(s(2))*(m1+m2) / l2*(m1+m2)-m2*l2*cos(s(1)-s(2))^2

This means that only the last component of tne numerator is divided, so a parenthesis is missing:

( ... - g*sin(s(2))*(m1+m2) ) / l2*(m1+m2)-m2*l2*cos(s(1)-s(2))^2
^                           ^

For the first term of the 3rd component you have

-m2 * l1 * (s(4))^2 ...

but in the publication there is:

-m2 * l1 * (s(3))^2 ...

Writing this formula in a a large unreadable line as an anonymous function makes it hard to debug. Don't do this. Prefer a dedicated function, which allows to split the parts and check them using the debugger:

function dF = DoublePendulum(t, s, m1, m2, l1, l2, g)
s12 = s(1) - s(2);
m12 = m1 + m2;

n1 = -m2 * l1 * s(3)^2 * sin(s12) * cos(s12) + g * m2 * sin(s(2)) * cos(s12) - ...
   m2 * l2 * s(4)^2 * sin(s12) - m12 * g * sin(s(1));
d1 = l1 * m12 - m2 * l1 * cos(s12)^2;

n2 =  m2 * l2 * s(4)^2 * sin(s12) * cos(s12) + g * m12 * sin(s(1)) * cos(s12) + ...
   m12 * l1 * s(3)^2 * sin(s12) - m12 * g * sin(s(2));
d2 = l2 * m12 - m2 * l2 * cos(s12)^2;

dF = [s(3); ...
   s(4); ...
   n1 / d1, ...
   n2 / d2];
end

Now the pattern of the two equations is visible on first sight and finding typos like the s(3) or s(4) is much easier.

Then define in your code:

F = @(t,s) DoublePendulum(t, s, m1, m2, l1, l2, g);

By the way: It is a drawback to obtain the dependent X/Y coordinates befor the integeration. Integrating the formulas for the angles is much faster and more accurate. Then the results can be used easily to compute the karthesian coordinates afterwards for the visualization. For a double pendulum both will work, but if you simulate > 4 connected pendulums the computing time with dependent coordinates explodes and the accuracy degrades.

[Jan]

No, for the matrix multiplication k*w, w must be a column vector. When your Runge Kutta method replies the same trajectories as Matlab's ODE45 as long as the step width is reduced, you can trust your integrator - and you should reduce the step width. Your RK is working fine.

Did you see the problem, that your ODE contained typos?

Answer 2

[John BG]

Hi Erlich

There's a double pendulum already implemented here:

https://uk.mathworks.com/matlabcentral/fileexchange/27212-animated-double-pendulum?s_tid=srchtitle

John BG

<mailto:[e-mail] [e-mail]>