结构动力学大作业

%

% plot(EWave.Time,EWave.Acel); %%

global MDOF %

MDOF.ND = 5; %

MDOF.MVec = zeros(MDOF.ND,1); MDOF.KVec = zeros(MDOF.ND,1);

for i = 1:5

MDOF.MVec(i) = 1000E3; end

for i = 1:5

MDOF.KVec(i) = 1.0*2000E6; end %%

MDOF.MMatrix = zeros(MDOF.ND, MDOF.ND); MDOF.KMatrix = zeros(MDOF.ND, MDOF.ND); %

for i = 1:MDOF.ND

MDOF.MMatrix(i,i) = MDOF.MVec(i); end %

for i = 1:MDOF.ND if i == 1

MDOF.KMatrix(i,i ) = MDOF.KVec(i ) + MDOF.KVec(i+1); MDOF.KMatrix(i,i+1) = -MDOF.KVec(i+1); else

if i == MDOF.ND

MDOF.KMatrix(i,i-1) = -MDOF.KVec(i); MDOF.KMatrix(i,i ) = MDOF.KVec(i); else

MDOF.KMatrix(i,i-1) = -MDOF.KVec(i );

MDOF.KMatrix(i,i ) = MDOF.KVec(i ) + MDOF.KVec(i+1); MDOF.KMatrix(i,i+1) = -MDOF.KVec(i+1); end end end %

[MDOF.Eigen_Vec, MDOF.Eigen_Val] = eig(MDOF.MMatrix\\MDOF.KMatrix); %

MDOF.GMMatrix=

transpose(MDOF.Eigen_Vec)*MDOF.MMatrix*MDOF.Eigen_Vec; MDOF.GKMatrix=

transpose(MDOF.Eigen_Vec)*MDOF.KMatrix*MDOF.Eigen_Vec; MDOF.GEVector=

transpose(MDOF.Eigen_Vec)*MDOF.MMatrix*ones(MDOF.ND,1); %%

global SDOF %

Result.Time=EWave.Time;

Result.Disp=zeros(length(Result.Time),MDOF.ND); Result.Vel=zeros(length(Result.Time),MDOF.ND); for i=1:MDOF.ND

SDOF.T=2*pirt(MDOF.Eigen_Val(i,i)); SDOF.W=2*pi/SDOF.T; SDOF.D=0.00; %

SDOF.AMatrix=zeros(2,2); %

SDOF.AMatrix(1,2)=1;

SDOF.AMatrix(2,1)=-2*SDOF.D*SDOF.W; SDOF.AMatrix(2,2)=-SDOF.W^2; %

SDOF.BVector=zeros(2,1); %

SDOF.BVector(2)=-MDOF.GEVector(i)/MDOF.GMMatrix(i,i); % Execute the time-history analysis

Solver.TD=[min(EWave.Time) max(EWave.Time)]; Solver.IC=zeros(2,1);

Solver.Opt=odeset('MaxStep',2.0*EWave.DT); %

[T,V]=ode45(@SDOF_Time_History_Analysis_ODE,Solver.TD,Solver.IC,Solver.Opt); %

for j=1:length(EWave.Time)

Result.Disp(j,i)=interp1(T,V(:,1),EWave.Time(j)); Result.Vel(j,i)=interp1(T,V(:,2),EWave.Time(j)); end %

disp(i); %

% subplot(2,1,1)

% plot(T, V(:,1),'blue','LineWidth',2); grid on; hold on; %

% subplot(2,1,2)

% plot(T, V(:,2),'red','LineWidth',2); grid on; hold on; end %%

Result.Disp_ALL=transpose(MDOF.Eigen_Vec*transpose( Result.Disp));

Result.Disp_1st=transpose(MDOF.Eigen_Vec(:,1)*transpose( Result.Disp(:,1)));

Result.Disp_2nd=transpose(MDOF.Eigen_Vec(:,1:2)*transpose( Result.Disp(:,1:2))); Result.Disp_3rd=transpose(MDOF.Eigen_Vec(:,1:3)*transpose( Result.Disp(:,1:3))); %

plot(Result.Time, Result.Disp_ALL(:,MDOF.ND),'red'); hold on;

plot(Result.Time, Result.Disp_1st(:,MDOF.ND),'blue'); hold on;grid on; 3、瑞利法程序ruilifa %% clc; %% clc; clear; close; %%

global MDOF %

MDOF.ND = 5; %

MDOF.MVec = zeros(MDOF.ND,1); MDOF.KVec = zeros(MDOF.ND,1);

for i = 1:5

MDOF.MVec(i) = 1000E3; end

for i = 1:5

MDOF.KVec(i) = 2000E6; end %%

MDOF.X = zeros(MDOF.ND,1); MDOF.X1 = zeros(MDOF.ND,1); for i=1:MDOF.ND s=0;

for n=i:MDOF.ND

s=s+MDOF.MVec(n); end

MDOF.X1(i)=s*9.8/(MDOF.KVec(i)); end

for i=1:MDOF.ND

if i==1

MDOF.X(i)=MDOF.X1(i); else

MDOF.X(i)=MDOF.X(i-1)+MDOF.X1(i); end end

MDOF.WVec=sqrt((9.8*transpose(MDOF.X)*MDOF.MVec)/(transpose(MDOF.MVec)*MDOF.X.^2));

4、矩阵迭代法Diedaifa %% clc; clear; close; %%

global MDOF %

MDOF.ND = 5; %

MDOF.MVec = zeros(MDOF.ND,1); MDOF.KVec = zeros(MDOF.ND,1);

for i = 1:5

MDOF.MVec(i) = 1000E3; end

for i = 1:5

MDOF.KVec(i) = 1.0*2000E6; end %%

MDOF.MMatrix = zeros(MDOF.ND, MDOF.ND); MDOF.KMatrix = zeros(MDOF.ND, MDOF.ND); %

for i = 1:MDOF.ND

MDOF.MMatrix(i,i) = MDOF.MVec(i); end %

for i = 1:MDOF.ND if i == 1

MDOF.KMatrix(i,i ) = MDOF.KVec(i ) + MDOF.KVec(i+1); MDOF.KMatrix(i,i+1) = -MDOF.KVec(i+1); else

if i == MDOF.ND

MDOF.KMatrix(i,i-1) = -MDOF.KVec(i); MDOF.KMatrix(i,i ) = MDOF.KVec(i);