clear all;
global M C K
% 输入质量矩阵，阻尼矩阵和刚度矩阵
a1=[153.0,170.0,170.0,170.0,170.0,170.0,183.0];
M=diag(a1);
b1=[3510,13500,14500,14600,14600,14600,13900];
b2=[-5780,-8950,-8740,-8720,-8710,-8430];
b3=[425,1090,1160,1170,1230];
b4=[1830,215,231,251];
b5=[121,69.4,78.5];
b6=[29.6,16.5];
b7=[8.74];
[C]=diag(b1)+diag(b2,1)+diag(b2,-1)+diag(b3,2)+diag(b3,-2)+diag(b4,3)+diag(b4,-3)+diag(b5,4)+diag(b5,-4)+diag(b6,5)+diag(b6,-5)+diag(b7,6)+diag(b7,-6);
c1=[204,808,870,874,875,876,833];
c2=[-350,-542,-529,-528,-527,-510];
c3=[111,66.2,70.3,70.7,74.6];
c4=[25.7,13.0,14.0,15.2,];
c5=[7.30,4.20,4.75];
c6=[1.79,1.00];
c7=[0.529];
[K]=diag(c1)+diag(c2,1)+diag(c2,-1)+diag(c3,2)+diag(c3,-2)+diag(c4,3)+diag(c4,-3)+diag(c5,4)+diag(c5,-4)+diag(c6,5)+diag(c6,-5)+diag(c7,6)+diag(c7,-6);
% newmark-b法
n=100;
t=0;
r=1/2;
b=1/4;
a0=1/(b*0.02^2);
a1=r/(b*0.02);
a2=1/(b*0.02);
a3=1/(2*b)-1;
a4=r/b-1;
a5=0.02/2*(r/b-2);
a6=0.02*(1-r);
a7=r*0.02;
u(:,1)=zeros(7,1);
u1(:,1)=zeros(7,1);
u2(:,1)=zeros(7,1);
p(:,1)=sin(pi*t)*[153.0,170.0,170.0,170.0,170.0,170.0,183.0]';
kk=[K]+a0*[M]+a1*[C];
for i=1:n
t=0.02*i;
p(:,i+1)=sin(pi*t)*[153.0,170.0,170.0,170.0,170.0,170.0,183.0]';
pp(:,i+1)=p(:,i+1)+[M]*(a0*u(:,i)+a2*u1(:,i)+a3*u2(:,i))+[C]*(a1*u(:,i)+a4*u1(:,i)+a5*u2(:,i));
u(:,i+1)=inv(kk)*pp(:,i+1);
u2(:,i+1)=-inv([M]+r*0.02*[C]+b*[K]*0.02*0.02)*(((1-r)*0.02*[C]+(1/2-b)*[K]*0.02*0.02)*u2(:,i)+([C]+0.02*[K])*u1(:,i)+[K]*u(:,i)-p(:,i+1));
u1(:,i+1)=u1(:,i)+(1-r)*0.02*u2(:,i)+r*0.02*u2(:,i+1);
u(:,i+1)=u(:,i)+0.02*u1(:,i)+(1/2-b)*0.02*0.02*u2(:,i)+b*0.02*0.02*u2(:,i+1);
end
X=[];
Y=[];
tt=[];
for i=1:n
tt=[tt,0.02*i];
X=[X,u(:,i)];
Y=[Y,tt(:,i)];
end
plot(t,u(:,i));