% Program 6.3 Animation program for pendulum using IVP solver
% Inputs: int = [a b] time interval,
%  initial values ic = [y(1,1) y(1,2)], h = step size
% Calls a one-step method such as trapstep.m
% Example usage: pend([0 10],[pi/2 0],.005)
function Damped_Pend_Dylan(int,ic,h)
n=round((int(2)-int(1))/h); % plot n points in total
y(1,:)=ic;                % enter initial conds in y
t(1)=int(1);
set(gca,'xlim',[-1.2 1.2],'ylim',[-1.2 1.2], ...
  'xtick',[-1 0 1],'ytick',[-1 0 1], ...
  'drawmode','fast','Visible','on')%,'NextPlot','add');
cla;
axis square               % make aspect ratio 1 - 1
bob = line('color','r','Marker','.','markersize',40,...
    'erase','xor','xdata',[],'ydata',[]);
rod = line('color','b','LineStyle','-','LineWidth',3,...
    'erase','xor','xdata',[],'ydata',[]);

MakeQTMovie('start','PendP4b.mov');
%filename is the name of the file
%we can do Q4M1.mov, Q4M2.mov, Q5M1.mov, or something or whatever you want
MakeQTMovie('framerate',110)
MakeQTMovie('quality',50);
h3 = text(-.5,.6,'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);
  xbob = cos(y(k+1,1)-pi/2);
  ybob = sin(y(k+1,1)-pi/2);
  abs(radtodeg(y(k+1,1)));

  % 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(-.5,.9,['Pendulum Angle: ',num2str(abs(radtodeg(y(k+1,1))))],...
  'BackgroundColor',[.7 .5 .7]);
  h2 = text(-.6,.75,['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)) < .01
      delete(h3);
      h3 = text(-.6,.6,...
          ['Angle at V = ZERO: ',num2str(abs(radtodeg(y(k+1,1))))],...
          'BackgroundColor',[.7 .9 .7]);
  end

  xrod = [0 xbob]; yrod = [0 ybob];
  set(rod,'xdata',xrod,'ydata',yrod)
  set(bob,'xdata',xbob,'ydata',ybob)
  drawnow; pause(h)

  % Delete's the text box objects displaying angle and velocity so
  % they don't overwrite on the figure and get all messy.


  %MakeQTMovie('addfigure');
  delete(h1);
  delete(h2);
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 = 1;%damping force
A = 10;
z(2) = -(g/length)*sin(y(1))-d*y(2);

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

Published with MATLAB® R2013a