Page 386 - MATLAB an introduction with applications

Direct Numerical Integration Methods ——— 371

M=[1 0;0 10];
C=[0.5 –0.1;–0.1 0.1];
dt=0.05;T=2;
X0=[0;0];X0d=[0;0];F0=[0;10];
X0dd=inv(M)*(F0–C*X0d–K*X0);
theta=1.4;
a0=6/(theta*dt)^2;a1=3/(theta*dt);a2=2*a1;
a3=2;a4=(1/2*theta*dt);a5=–a2/theta;
a6=1–3/theta;
a7=dt/2;a8=dt^2/6;
Kb=K+a0*M+a1*C;
i=1;
X(:,1)=X0;Xd(:,1)=X0d;Xdd(:,1)=X0dd;t=0;
fprintf(‘time(s)\t\tX1\t\tX2\n’);
fprintf(‘%f\t%f\t%f\n’,t,X(1,1),X(2,1));
for t=dt:dt:T
i=i+1;
F=[0;10];
Ftb=F0+M*(a0*X(:,i–1)+a2*Xd(:,i–1)+a3*Xdd(:,i–1))+C*(a1*X(:,i–1)
+ a3*Xd(:,i–1)+a4*Xdd(:,i–1))+theta*(F–F0);
Xt(:,i)=inv(Kb)*Ftb;
Xdd(:,i)=(a0/theta)*(Xt(:,i)–X(:,i–1))+a5*Xd(:,i–1)+a6*Xdd(:,i–1);
Xd(:,i)=Xd(:,i–1)+a7*(Xdd(:,i)+Xdd(:,i–1));
X(:,i)=X(:,i–1)+dt*Xd(:,i–1)+a8*(Xdd(:,i)+2*Xdd(:,i–1));
F0=F;
fprintf(‘%f\t%f\t%f\n’,t,X(1,i),X(2,i));
end
t=[0:dt:T];
plot(t,X(1,:),’–p’,t,X(2,:),’–*’)
xlabel(‘time(s)’);
ylabel(‘displacement(m)’);
legend(‘DOF–1’,’DOF–2');
grid on;

The output is shown below:
time X(1) X(2)
0.000000 0.000000 0.000000
0.050000 0.000003 0.001250
0.100000 0.000022 0.004998