Page 384 - MATLAB an introduction with applications

Direct Numerical Integration Methods ——— 369

X(:,1)=X0;
Fbar=F–kbar*X0–cbar*Xprev;
X(:,2)=inv(mbar)*Fbar;
Fbar=F–kbar*X(:,2)–cbar*X0;
X(:,3)=inv(mbar)*Fbar;

% HOUBOLT METHOD BEGINS
a0=2/(dt^2);a1=11/(6*dt);a2=5/(dt^2);a3=3/dt;a4=–2*a0;
a5=–a3/2;a6=a0/2;a7=a3/9;
Kb=K+a0*M+a1*C;
p=3;
for i=3:length(t)
F=[0;10];% F(t+2dt)
Fb=F+M*(a2*X(:,i)+a4*X(:,i–1)+a6*X(:,i–2))+ C*(a3*X(:,i)+a5*X(:,i–
1)+a7*X(:,i–2));
X(:,i+1)=inv(Kb)*Fb;
Xdd(:,i+1)=a0*X(:,i+1)–a2*X(:,i)–a4*X(:,i–1)–a6*X(:,i–2);
Xd(:,i+1)=a1*X(:,i+1)–a3*X(:,i)–a5*X(:,i–1)–a7*X(:,i–2);
p=p+1;
end
fprintf(‘\ntime\t\tX1\t\tX2\n’);
for i=1:p
time(i)=(i–1)*dt;
fprintf(‘%f\t%f\t%f\n’,time(i),X(1,i),X(2,i))
end
plot(time,X(1,:), ‘–p’,time,X(2,:),‘–*’);
grid on;
xlabel(‘time(s)’);
ylabel(‘displacement (m)’);
legend(‘DOF-1’,’DOF-2');

The output of the program is given below:
time X(1) X(2)
0.000000 0.000000 0.000000
0.050000 –0.000000 0.001250
0.100000 0.000015 0.004998
0.150000 0.000069 0.011243
0.200000 0.000184 0.019980
0.250000 0.000387 0.031207
0.300000 0.000704 0.044920