可是我用newmark法求出来的响应图片怎么是这样的啊
K=sysK;% K ----- 调用刚度矩阵
M=sysM; % M ----- 调用质量矩阵
%C=0.03*M+0.02*K;% C ----- 阻尼矩阵
C=0;
dt=0.01;% dt ----- 时间步长
tend=20;% tend --- 结束时间
w=600;
[n,n] = size( K ) ;
%Newmark算法
gama = 0.5 ;
beta = 0.25 ;
[n,n] = size( K ) ;
Nalpha0 = 1/beta/dt^2 ;
Nalpha1 = gama/beta/dt ;
Nalpha2 = 1/beta/dt ;
Nalpha3 = 1/2/beta - 1 ;
Nalpha4 = gama/beta - 1 ;
Nalpha5 = dt/2*(gama/beta-2) ;
Nalpha6 = dt*(1-gama) ;
Nalpha7 = gama*dt ;
NK1 = K + Nalpha0*M + Nalpha1*C ;
Nd = zeros( n, floor(tend/dt) + 1 ) ;
Nv = zeros( n, floor(tend/dt) + 1 ) ;
Na = zeros( n, floor(tend/dt) + 1 ) ;
Nd(:,1) = 0.01 ; %%%%%%%%%%%%%%%%%%%%%%%%初始位移
Nv(:,1) = 0 ; %%%%%%%%%%%%%%%%%%%%%%%%初始速度
f(:,1) =0 ;%初始载荷
Na(:,1) = M\(f(:,1)-K*Nd(:,1)-C*Nv(:,1)); %初始加速度
t=0:dt:tend;
for i=2:1:length(t)
f(:,i)=0*100*sin(w*t(i));
f2 = f(:,i) + M*(Nalpha0*Nd(:,i-1)+Nalpha2*Nv(:,i-1)+Nalpha3*Na(:,i-1))+ C*(Nalpha1*Nd(:,i-1)+Nalpha4*Nv(:,i-1)+Nalpha5*Na(:,i-1)) ;
Nd(:,i) = inv(NK1)*f2 ;
Na(:,i) = Nalpha0*(Nd(:,i)-Nd(:,i-1)) - Nalpha2*Nv(:,i-1) - Nalpha3*Na(:,i-1) ;
Nv(:,i) = Nv(:,i-1) + Nalpha6*Na(:,i-1) + Nalpha7*Na(:,i) ;
end
figure ;
plot(t,Nd(1,:),'g')
xlabel('t');
ylabel('u(t)');
title('Newmark-beta位移与时间的关系');
|