%用Wilson方法計算結構動力反應(無阻尼)
clear;
m=[61.7 0;
   0 29.9];
k=[34016 -20092;
   -20092 20092];
xita=1.4;     %確定系數xita的值
n=length(m);
tt=3;             
pbc=zeros(n-1,1);
p(:,1)=[pbc;115.2];    %初始荷載、位移、速度、加速度
y(:,1)=zeros(n,1);
y1(:,1)=zeros(n,1);
y2(:,1)=inv(m)*(p(:,1)-k*y(:,1));
detat=0.02;            %計算時間步長
a0=6/((xita*detat)^2);
a1=3/(xita*detat);
a2=2*a1;
a3=xita*detat/2;
a4=a0/xita;
a5=-a2/xita;
a6=1-3/xita;
a7=detat/2;
a8=detat^2/6;
k2=k+a0*m;
for t=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
   r(:,j)=p(:,j-1)+xita*(p(:,j)-p(:,j-1))+m*(a0*y(:,j-1)+a2*y1(:,j-1)+2*y2(:,j-1));
   yy(:,j)=inv(k2)*r(:,j);
   y2(:,j)=a4*(yy(:,j)-y(:,j-1))+a5*y1(:,j-1)+a6*y2(:,j-1);
   y1(:,j)=y1(:,j-1)+a7*(y2(:,j)+y2(:,j-1));
   y(:,j)=y(:,j-1)+detat*y1(:,j-1)+a8*(y2(:,j)+2*y2(:,j-1));
   y2(:,j)=inv(m)*(p(:,j)-k*y(:,j));
end
tt2=0:detat:tt;
subplot(2,2,1),plot(tt2,y(1,:))
title('The Wilson Method')
xlabel('t')
ylabel('X1')
subplot(2,2,2),plot(tt2,y(2,:))
title('The Wilson Method')
xlabel('t')
ylabel('X2')
fid=fopen('shuju.m','w+');
fprintf(fid,'\nt       x1        x2\n');
for i=1:j
   fprintf(fid,'\n% +7.5f   %+12.4e   %12.4e',tt2(i),y(1,i),y(2,i));
end
fclose(fid);