How to solve the varying length of strings of double pendulum in the simulation?

In my code,while the simulation is running,the length of the strings vary though the lengths are constant in the code. How can I solve it?

 function double_pendulum(m1,m2,l1,l2,theta1,theta2,tmax)
%value of the constants
  m1=0.1;
  m2=0.1;
  l1=0.5;
  l2=0.5;
  theta1=86;
  theta2=86;
  g=9.8;
  tmax=20;
  t(1)=0;
 %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(3))^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;

   while t(i)<tmax %==================================================
 %Computation of position and angular Velocity
   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,:));
    arraysize=size(t);
    axissize = [min(min(x1),min(x2)) max(max(x1),max(x2)) min(min(y1),min(y2)) max(max(y1),max(y2))];  % Sets axis size for animation.
    max(arraysize);
    pendulumtopx = sum(x1)/max(arraysize);
    pendulumtopy = sum(y1)/max(arraysize);

   for i = 1 : max(arraysize)

    plotarrayy = [pendulumtopy y1(i)];  % Plots solution at each time interval
    plotarrayx = [pendulumtopx x1(i)];
    plotarray2x = [x1(i) x2(i)];
    plotarray2y = [y1(i) y2(i)];
    plot(x1(i),y1(i),'o',x2(i),y2(i),'o','markersize',10,'markerfacecolor','g')
    hold on
    plot(x1(1:i),y1(1:i),'r');
    plot(x2(1:i),y2(1:i),'b');
    plot(plotarrayx,plotarrayy)
    plot(plotarray2x,plotarray2y)       
    title(['Double pendulum simulation']) 
    hold off

   xlabel('x','fontsize',12)   
   ylabel('y','fontsize',12)
   axis(axissize);
    pause(0.001);   % Pause time between animations.
   i=i+1;
end

 figure
 subplot(4,1,1)
 hold on
 plot(s(1,:),s(2,:),'b')

 title('Double pendulum phase portrait','fontsize',12)
        xlabel('theta1','fontsize',12)
        ylabel('theta2','fontsize',12)
 subplot(4,1,2)
 hold on
 plot(s(1,:),s(3,:),'r')
 title('Double pendulum phase portrait','fontsize',12)
        xlabel('theta1','fontsize',12)
        ylabel('velocity1','fontsize',12)
 %        axis([min(s1) max(s1) min(s2) max(s2)])

 subplot(4,1,3)% Plotting angle velocity vs time
 plot(t,s(3,:),t,0)
 title('Angle velocity vs Time Graph');
 xlabel('time');
 ylabel('\theta1(t)')
 subplot(4,1,4)
 plot(t,s(4,:),t,0)
 set(gca,'XLim',[0 tmax])
 title('Angle velocity vs Time Graph')
 xlabel('time')
 ylabel('\theta2(t)')

Answer 1

[James Tursa]

Do you have a profile of how the string lengths change as a function of time? E.g., can you do something like this:

Have functions defined

function L = l1(t)
L = % put code here for length as function of t
end

function L = l2(t)
L = % put code here for length as function of t
end

And then make your function handle call l1 and l2 as functions of t:

F=@(t,s)[s(3) s(4) (-m2*l1(t)*(s(3))^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(t)*(m1+m2)-m2*l1(t)*(cos(s(1)-s(2))^2) (m2*l2(t)*(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(t)*(s(3))^2*sin(s(1)-s(2))*(m1+m2)-g*sin(s(2))*(m1+m2))/l2(t)*(m1+m2)-m2*l2(t)*cos(s(1)-s(2))^2]';

Of if the lengths are a function of something else, put that into the l1 and l2 function code.

[reema shrestha]

I tried changing the value of axissize and I think it worked. May be the problem arose because of it.