Page 312 - Applied Numerical Methods Using MATLAB

PROBLEMS 301

%do_MBK
clf
t0=0;tf=10; x0=[00];
[t1,x] = ode_Ham(’f_MBK’,[t0 tf],x0);
dt = t1(2) - t1(1);
for n = 1:length(t1)
u(n) = udu_MBK(t1(n));
end
figure(1), clf
animation = 1;
if animation
figure(2), clf
draw_MBK(5,1,x(1,2),u(1))
axis([-2 2 -1 14]), axis(’equal’)
pause
for n = 1:length(t1)
clf, draw_MBK(5,1,x(n,2),u(n),’b’)
axis([-2 2 -1 14]), axis(’equal’)
pause(dt)
figure(1)
plot(t1(n),u(n),’r.’, t1(n),x(n,2),’b.’)
axis([0 tf -0.2 1.2]), hold on
figure(2)
end
draw_MBK(5,1,x(n,2),u(n))
axis([-2 2 -1 14]), axis(’equal’)
end
function [u,du] = udu_MBK(t)
i = fix(t);
if mod(i,2) == 0, u = t-i; du = 1;
elseu=1- t+i;du=-1;
end
function draw_MBK(n,w,y,u,color)
%n: the # of spring windings
%w: the width of each object
%y: displacement of the top of MBK
%u: displacement of the bottom of MBK
if nargin < 5, color = ’k’; end
p1 = [-w u + 4]; p2 = [-w 9 + y];
xm = 0; ym = (p1(2) + p2(2))/2;
xM = xm + w*1.2*[-1 -1 1 1 -1];
yM = p2(2) + w*[1 3 3 1 1];
plot(xM,yM,color), hold on %Mass
spring(n,p1,p2,w,color) %Spring
damper(xm + w,p1(2),p2(2),w,color) %Damper
wheel_my(xm,p1(2)- 3*w,w,color) %Wheel

function dx = f_MBK(t,x)
M = 1; B = 0.1; K = 0.1;
[u,du] = udu_MBK(t);
dx = x*[0 1; -B/M - K/M]’+[0 (K*u + B*du)/M];