tacoma

 Program 6.6 Animation program for bridge using IVP solver
 Inputs: inter = [a b] time interval, 
   ic = [y(1,1) y(1,2) y(1,3) y(1,4)],
   h = stepsize, p = steps per point plotted
 Calls a one-step method such as trapstep.m
 Example usage: tacoma([0 1000],[0 0 0.001 0],.04,3);
 clf                                  % clear figure window

0001 % Program 6.6 Animation program for bridge using IVP solver
0002 % Inputs: inter = [a b] time interval,
0003 %   ic = [y(1,1) y(1,2) y(1,3) y(1,4)],
0004 %   h = stepsize, p = steps per point plotted
0005 % Calls a one-step method such as trapstep.m
0006 % Example usage: tacoma([0 1000],[0 0 0.001 0],.04,3);
0007 function [height, theta] = tacoma(inter,ic,h,p)
0008 % clf                                  % clear figure window
0009 a=inter(1);b=inter(2);n=ceil((b-a)/(h*p));% plot n points
0010 % y(1) = y; y(2) = y'; y(3) = theta; y(4) = theta'
0011 y(1,:)=ic;                           % enter initial conds in y
0012 t(1)=a;len=6;
0013 % set(gca,'XLim',[-8 8],'YLim',[-8 8], ...
0014 %    'XTick',[-8 0 8],'YTick',[-8 0 8], ...
0015 %    'Drawmode','fast','Visible','on','NextPlot','add');
0016 % cla;                                 % clear screen
0017 % axis square                          % make aspect ratio 1 - 1
0018 % road=line('color','b','LineStyle','-','LineWidth',5,...
0019 %     'erase','xor','xdata',[],'ydata',[]);
0020 % lcable=line('color','r','LineStyle','-','LineWidth',1,...
0021 %     'erase','xor','xdata',[],'ydata',[]);
0022 % rcable=line('color','r','LineStyle','-','LineWidth',1,...
0023 %     'erase','xor','xdata',[],'ydata',[]);
0024 % pause on;
0025 height = zeros(n*p,1);
0026 theta = zeros(n*p,1);
0027 for k=0:n-1 
0028   for i=1:p
0029     t(i+1) = t(i)+h;
0030 %     y(i+1,:) = trapstep(t(i),y(i,:),h);
0031     y(i+1,:) = rk4step(t(i),y(i,:),h);
0032   end
0033   
0034   height(k*p+1:k*p+p+1) = y(:,1);
0035   theta(k*p+1:k*p+p+1) = y(:,3);
0036   
0037   y(1,:) = y(p+1,:);
0038   t(1)=t(p+1);
0039   z1(k+1)=y(1,1);
0040   z3(k+1)=y(1,3);
0041   c=len*cos(y(1,3));s=len*sin(y(1,3));
0042 %   set(road,'xdata',[-c c],'ydata',[-s-y(1,1) s-y(1,1)])
0043 %   set(lcable,'xdata',[-c -c],'ydata',[-s-y(1,1) 8])
0044 %   set(rcable,'xdata',[c c],'ydata',[s-y(1,1) 8])
0045 %   drawnow;
0046 %   pause(h);
0047 end
0048 pause off;
0049 figure(100);
0050 plot(height);
0051 figure(101);
0052 plot(theta);
0053 
0054 function y = trapstep(t,x,h)
0055 %one step of the Trapezoid Method
0056 z1=ydot(t,x);
0057 g=x+h*z1;
0058 z2=ydot(t+h,g);
0059 y=x+h*(z1+z2)/2;
0060 
0061 function y=rk4step(t,w,h)
0062 %one step of the Runge-Kutta order 4 method
0063 s1=ydot(t,w);
0064 s2=ydot(t+h/2,w+h*s1/2);
0065 s3=ydot(t+h/2,w+h*s2/2);
0066 s4=ydot(t+h,w+h*s3);
0067 y=w+h*(s1+2*s2+2*s3+s4)/6;
0068 
0069 function ydot=ydot(t,y)
0070 len=6; a=0.2; W=80; omega=2*pi*38/60;
0071 a1=exp(a*(y(1)-len*sin(y(3))));
0072 a2=exp(a*(y(1)+len*sin(y(3)))); 
0073 ydot(1) = y(2);
0074 ydot(2) = -0.01*y(2)-0.4*(a1+a2-2)/a+0.2*W*sin(omega*t);
0075 ydot(3) = y(4);
0076 ydot(4) = -0.01*y(4)+1.2*cos(y(3))*(a1-a2)/(len*a);