Page 362 - MATLAB an introduction with applications

Direct Numerical Integration Methods ——— 347

F0=0;F=150;
T=1;
x0dd=inv(m)*(F0-c*x0d-k*x0);
xprev=x0-(dt*x0d)+((dt^2)*x0dd/2);
a0=1/dt^2;a1=1/(2*dt);a2=2*a0;
mbar=a0*m+a1*c;
t=0;
v(1)=x0d;a(1)=x0dd;
i=1;
for t=0:dt:T+dt
X(i)=x0;
if t<=0.2 f=(F*t/0.2);
else if (t>0.2 & t<=0.4) f=-(F/0.2).*(t-0.4);
else if t>0.4 f=0;
end
end
end
fbar=f+(a2*m-k)*x0+(a1*c-a0*m)*xprev;
x=inv(mbar)*fbar;
xprev=x0;
x0=x;
i=i+1;
p=i;
end
for i=2:p-1
if i<p-1
v(i)=(X(i+1)-X(i-1))/(2*dt);
a(i)=(X(i+1)-2*X(i)+X(i-1))/dt^2;
end
end
fprintf(‘\ntime\t\tdisplacement\tvelocity\tacceleration\n’);
i=1;
for t=0:dt:T
fprintf(‘%f\t%f\t%f\t%f\n’,t,X(i),v(i),a(i));
i=i+1;
end
t=[0:dt:T+dt];
plot(t,X,’-p’);
xlabel(‘time(s)’);
ylabel(‘displacement(m)’);
grid on;