%用中心差分法計算結構動力反應(無阻尼)
m=[61.7 0;
   0 29.9];
k=[34016 -20092;
   -20092 20092];
n=length(m);
tt=3;             
pbc=zeros(n-1,1);
p(:,1)=[pbc;115.2];%初始荷載、位移、速度、加速度
p(:,2)=[pbc;115.2];
z(:,1)=zeros(n,1);
z1(:,1)=zeros(n,1);
z2(:,1)=inv(m)*(p(:,1)-k*z(:,1));
detat=0.02;            %計算時間步長
a0=1/(detat^2);
a1=1/(2*detat);
zz(:,1)=z(:,1)-z1(:,1)/(2*a1)+z2(:,1)/(2*a0);
z(:,2)=inv(a0*m)*(p(:,1)-(k-2*a0*m)*z(:,1)-a0*m*zz(:,1));
for t=2*detat:detat:tt
   j=round(t/detat)+1;
   if t<=0.4 
      p(:,j)=[pbc;115.2];   %定義變化的激勵
   elseif t>0.4
      p(:,j)=[pbc;0];
   end
   z(:,j)=inv(a0*m)*(p(:,j-1)-(k-2*a0*m)*z(:,j-1)-a0*m*z(:,j-2));
end
tt2=0:detat:tt;
subplot(2,2,1),plot(tt2,z(1,:))
title('中心差分法')
xlabel('t')
ylabel('X1')
subplot(2,2,2),plot(tt2,z(2,:))
title('中心差分法')
xlabel('t')
ylabel('X2')
fid=fopen('shuju2.m','w+');
fprintf(fid,'\nt       x1        x2\n');
for i=1:j
   fprintf(fid,'\n% +7.5f   %+12.4e   %12.4e',tt2(i),z(1,i),z(2,i));
end
fclose(fid);