%********************************************************************%子空間迭代%********************************************************************FP1=fopen('input.txt','rt');%打開數據輸入文本niter=fscanf(FP1,'%f',1);%讀入最多迭代次數nitern=fscanf(FP1,'%f',1);%讀入自由度數目ncolumn=fscanf(FP1,'%f',1);%讀入假設振型階數columnerr=fscanf(FP1,'%f',1);%讀入誤差限errA0=xlsread('假設振型');%從excle讀入假設振型M=xlsread('質量矩陣');%從excle讀入質量矩陣K=xlsread('剛度矩陣');%從excle讀入剛度矩陣forj=1:columnA0(:,j)=A0(:,j)/max(abs(A0(:,j)));%將假設振型A0基準化end%----開始迭代foriter=1:niter%----第iter次迭代F0=A0;F1=K\M*F0;forj=1:columnF1(:,j)=F1(:,j)/max(abs(F1(:,j)));endM1=F1'*M*F1;K1=F1'*K*F1;[V1,W1]=eig(M1\K1);forj=1:columna1(:,j)=V1(:,j)/V1(j,j);endA1=F1*a1;forj=1:columnA1(:,j)=A1(:,j)/max(abs(A1(:,j)));end%----精度判斷fori=1:nforj=1:columnifabs(A1(i,j)-A0(i,j))>=errA0=A1;elseifabs(A1(i,j)-A0(i,j))<errbreakendendendendendW=sqrt(diag(W1));A=A1;%----輸出自振特性xlswrite('自振特性.xlsx',A,1,'A2');%將數據寫入xlsx文件中，數據的開始位置為sheet1，A2xlswrite('自振特性.xlsx',W,2,'A2');%將數據寫入xlsx文件中，數據的開始位置為sheet2，A2%**************************************************************************%振型疊加法%**************************************************************************FP1=fopen('input.txt','rt');%打開數據輸入文本dt=fscanf(FP1,'%f',1);%讀入時間步長dtnt=fscanf(FP1,'%f',1);%讀入計算次數ntn=fscanf(FP1,'%f',1);%讀入自由度數目nnmode=fscanf(FP1,'%f',1);%讀入選取振型數目nmodep=fscanf(FP1,'%f',1);%讀入外力荷載參數wbar=fscanf(FP1,'%f',1);%讀入外力荷載頻率rdamp=xlsread('阻尼比');%從excle讀入阻尼比M=xlsread('質量矩陣');%從excle讀入質量矩陣K=xlsread(剛度矩陣');%從excle讀入剛度矩陣A=xlsread('振型');%從excle讀入前nmode階振型W=xlsread('頻率');%從excle讀入前nmode階頻率%----確定與正則坐標對應的廣義荷載向量fori=1:nt;forj=1:n;ifi>=0&i<=nt;Q(j,i)=p*sin(i*dt*wbar);end;end;end;%----將振型進行正則化處理Mgnrl=diag(A'*M*A);fori=1:nmodeforj=1:nAbar(j,i)=A(j,i)/sqrt(Mgnrl(i));endend%----本程序以例4-5為實例，基于其正則坐標響應的解析式，寫出其穩態響應的程序，對于無解析式的荷載工況，此處需改為杜哈美數值積分程序fori=1:nmodeD(i)=1/sqrt((1-wbar^2/W(i)^2)^2+(2*rdamp(i)*wbar/W(i))^2);cta(i)=atan((2*rdamp(i)*wbar/W(i))/(1-wbar^2/W(i)^2));endfori=1:nmodeforj=1:ntTbar(i,j)=p*Abar(:,i)'*ones(n,1)*D(i)*sin(j*wbar*dt-cta(i))/W(i)^2;endend%----體系的穩態響應fori=1:nforj=1:ntd(i,j)=Abar(i,:)*Tbar(:,j);endend%----輸出穩態響應i=1:nt;xlswrite('穩態響應.xlsx',d(:,i),1,'A1');%將數據寫入xlsx文件中，數據的開始位置為sheet1，A1%**************************************************************************%newmark-β法%**************************************************************************FP1=fopen('input.txt','rt');%打開數據輸入文件dt=fscanf(FP1,'%f',1);%讀入時間步長dtnt=fscanf(FP1,'%f',1);%讀入計算次數ntn=fscanf(FP1,'%f',1);%讀入自由度數目nalfa=fscanf(FP1,'%f',1);%讀入積分控制參數alfadta=fscanf(FP1,'%f',1);%讀入積分控制參數dtaRayleigha0=fscanf(FP1,'%f',1);%讀入瑞利阻尼待定常數Rayleigha0Rayleigha1=fscanf(FP1,'%f',1);%讀入瑞利阻尼待定常數Rayleigha1p=fscanf(FP1,'%f',1);%讀入外力荷載參數wbar=fscanf(FP1,'%f',1);%讀入外力荷載頻率M=xlsread('質量矩陣');%從excle上讀入質量矩陣K=xlsread('剛度矩陣');%從excle上讀入剛度矩陣C=Rayleigha0*M+Rayleigha1*K;%定義阻尼矩陣C%----外力荷載fori=1:ntforj=1:nifi>0&&i<=ntQ(j,i)=p*sin(i*dt*wbar);endendend%----計算動力響應b0=1/(alfa*dt^2);b1=dta/(alfa*dt);b2=1/(alfa*dt);b3=1/(2*alfa)-1;b4=dta/alfa-1;b5=dt/2*(dta/alfa-2);d(:,1)=zeros(n,1);v(:,1)=zeros(n,1);a(:,1)=M\(Q(:,1)-C*v(:,1)+K*d(:,1));Keff=b0*M+b1*C+K;fori=2:ntQeff=Q(:,i)+M*(b0*d(:,i-1)+b2*v(:,i-1)+b3*a(:,i-1))+C*(b1*d(:,i-1)+b4*v(:,i-1)+b5*a(:,i-1));d(:,i)=Keff\Qeff;a(:,i)=b0*(d(:,i)-d(:,i-1))-b2*v(:,i-1)-b3*a(:,i-1);v(:,i)=b1*(d(:,i)-d(:,i-1))-b4*v(:,i-1)-b5*a(:,i-1);end%----輸出動力響應i=1:nt;xlswrite('動力響應.xlsx',d(:,i),1,'A1');%將數據寫入xlsx文件中，數據的開始位置為sheet1，A1xlswrite('動力響應.xlsx',v(:,i),2,'A1');%將數據寫入xlsx文件中，數據的開始位置為sheet2，A1xlswrite('動力響應.xlsx',a(:,i),3,'A1');%將數據寫入xlsx文件中，數據的開始位置為sheet3，A1%**************************************************************************%wilson-θ法%**************************************************************************FP1=fopen('input.txt','rt');%打開數據輸入文件dt=fscanf(FP1,'%f',1);%讀入時間步長dtnt=fscanf(FP1,'%f',1);%讀入計算次數ntn=fscanf(FP1,'%f',1);%讀入自由度數目ncta=fscanf(FP1,'%f',1);%讀入參數ctaRayleigha0=fscanf(FP1,'%f',1);%讀入瑞利阻尼待定常數Rayleigha0Rayleigha1=fscanf(FP1,'%f',1);%讀入瑞利阻尼待定常數Rayleigha1p=fscanf(FP1,'%f',1);%讀入外力荷載參數wbar=fscanf(FP1,'%f',1);%讀入外力荷載頻率M=xlsread('質量矩陣');%從excle上讀入質量矩陣K=xlsread('剛度矩陣');%從excle上讀入剛度矩陣C=Rayleigha0*M+Rayleigha1*K;%定義阻尼矩陣C%----外力荷載fori=1:ntforj=1:nifi>0&&i<=ntQ(j,i)=p*sin(i*dt*wbar);endendend%----計算動力響應b0=6/(cta^2*dt^2);b1=3/(cta*dt);b2=6/(cta*dt);b3=cta*dt/2;b4=6/(cta^3*dt^2);b5=-6/(cta^2*dt);b6=1-3/cta;b7=dt/2;b8=dt^2/6;d(:,1)=zeros(n,1);v(:,