% Program 6.3 Animation program for pendulum
% Inputs: time interval inter,
%  initial values ic = [y(1,1) y(1,2)], number of steps n
% Calls a one-step method such as trapstep.m
% Example usage: pend([0 10],[pi/2 0],200)
function Double_Pend_Don(inter,ic,n)
h=(inter(2)-inter(1))/n;  % plot n points in total
y(1,:)=ic;                % enter initial conds in y
t(1)=inter(1);
%xxx = zeros(10);
%yyy = zeros(10);

subplot(1,2,1);

set(gca,'xlim',[-2.4 2.4],'ylim',[-2.4 2.4], ...
  'XTick',[-1 0 1],'YTick',[-1 0 1], ...
  'Drawmode','fast','Visible','on','NextPlot','add');
cla;
axis square

subplot(1,2,2);

set(gca,'xlim',[-5 5],'ylim',[-5 5], ...
  'XTick',[-1 0 1],'YTick',[-1 0 1], ...
  'Drawmode','fast','Visible','on','NextPlot','add');
cla;
axis square

scatter([0],[0]);

subplot(1,2,1);


bob1 = line('color','r','Marker','.','markersize',40,...
    'erase','xor','xdata',[],'ydata',[]);
rod1 = line('color','b','LineStyle','-','LineWidth',3,...
    'erase','xor','xdata',[],'ydata',[]);
bob2 = line('color','r','Marker','.','markersize',40,...
    'erase','xor','xdata',[],'ydata',[]);
rod2 = line('color','b','LineStyle','-','LineWidth',3,...
    'erase','xor','xdata',[],'ydata',[]);

MakeQTMovie('start','filename.mov'); %filename is the name of the file
                                     %we can do Q4M1.mov, Q4M2.mov, Q5M1.mov, or something
MakeQTMovie('framerate',20); %or whatever you want

fig = figure(1);
set(fig, 'Position', [100, 100, 1100, 1100*9/16]);

h3 = text(-2.2,1.1,'Angle at V = ZERO: N/A',...
'BackgroundColor',[.7 .9 .7]);

h6 = text(.3,1.1,'Angle at V = ZERO: N/A',...
'BackgroundColor',[.7 .9 .7]);

for k=1:n


  t(k+1)=t(k)+h;
  y(k+1,:)=trapstep(t(k),y(k,:),h);
  xbob1 = sin(y(k+1,1)); ybob1 = -cos(y(k+1,1));
  xrod1 = [0 xbob1]; yrod1 = [0 ybob1];
  xbob2 = xbob1 + sin(y(k+1,3)); ybob2 = ybob1 + -cos(y(k+1,3));
  xrod2 = [xbob1 xbob2]; yrod2 = [ybob1 ybob2];
  set(rod1,'xdata',xrod1,'ydata',yrod1)
  set(bob1,'xdata',xbob1,'ydata',ybob1)
  set(rod2,'xdata',xrod2,'ydata',yrod2)
  set(bob2,'xdata',xbob2,'ydata',ybob2)

  %plot(xbob2,ybob2);


    % Using the title to display the angle autocenters the text,
  % making it move around a lot, which didn't look great.
  %
  % The following code creates objects for the text boxes and manipulate's
  % the objects' properties.

  h1 = text(-2.2,1.65,['Pendulum Angle: ',num2str(abs(radtodeg(y(k+1,1))))],...
  'BackgroundColor',[.7 .5 .7]);
  h2 = text(-2.3,1.375,['Pendulum Velocity: ',num2str(abs(radtodeg(y(k+1,2))))],...
  'BackgroundColor',[.3 .8 .8]);

  % This if statement displays the max angle reached each time the
  % pendulum's velocity hits zero.
  if abs(y(k+1,2)) < .1
      delete(h3);
      h3 = text(-2.2,1.1,...
          ['Angle at V = ZERO: ',num2str(abs(radtodeg(y(k+1,1))))],...
          'BackgroundColor',[.7 .9 .7]);
  end

  %['Angle at V = ZERO: ',num2str(abs(radtodeg(y(k+1,1))))],...

  %radtodeg(y(k+1,1))

  % And this following code helps do it for the double pendulum! Yay!
    h4 = text(.3,1.65,['Pendulum Angle: ',num2str(abs(radtodeg(y(k+1,3))))],...
  'BackgroundColor',[.7 .5 .7]);
  h5 = text(.2,1.375,['Pendulum Velocity: ',num2str(abs(radtodeg(y(k+1,4))))],...
  'BackgroundColor',[.3 .8 .8]);

  % This if statement displays the max angle reached each time the
  % pendulum's velocity hits zero. Kinda like before.
  if abs(y(k+1,4)) < .1
      delete(h6);
      h6 = text(.3,1.1,...
          ['Angle at V = ZERO: ',num2str(abs(radtodeg(y(k+1,3))))],...
          'BackgroundColor',[.7 .9 .7]);
  end

  subplot(1,2,2);
    hold off;
    hold on;
 xxx(k) = .5*abs(y(k+1,2))*sin(y(k+1,1)) ;
 yyy(k) = -.5*abs(y(k+1,2))*cos(y(k+1,1));
 plot(xxx,yyy,'g');

 xxxx(k) = .5*abs(y(k+1,4))*sin(y(k+1,3)) ;
 yyyy(k) = -.5*abs(y(k+1,4))*cos(y(k+1,3));
 plot(xxxx,yyyy,'r');

 subplot(1,2,1);

  drawnow; pause(h)
  MakeQTMovie('addfigure');
  delete(h1);
  delete(h2);
  delete(h4);
  delete(h5);
end

MakeQTMovie('finish');

function y = trapstep(t,x,h)
%one step of the Trapezoid Method
z1=ydot(t,x);
g=x+h*z1;
z2=ydot(t+h,g);
y=x+h*(z1+z2)/2;

function z=ydot(t,y)
g=9.81;length=1;
%z(1) = y(2);
%z(2) = -(g/length)*sin(y(1));
d = .5;
z(1) = y(2);
z(2) = (-3*g*sin(y(1)) - g*sin(y(1) - 2*y(3)) - 2*sin(y(1) - y(3))*(y(4)^2 - (y(2)^2)*cos(y(1) - y(3))))  /  (3 - cos(2*y(1) - 2*y(3)))  -  d*y(2);
z(3) = y(4);
z(4) = (2*sin(y(1) - y(3))*(2*y(2)^2 + 2*g*cos(y(1)) + (y(4)^2)*cos(y(1)-y(3)))) / (3 - cos(2*y(1) - 2*y(3)));

Error using Double_Pend_Don (line 7)
Not enough input arguments.

Published with MATLAB® R2013a