newmark2.m

利用newmark方法计算无阻尼系统

clear;
%用Newmark方法计算结构动力反应(无阻尼)
m=[61.7 0;
   0 29.9];
k=[34016 -20092;
   -20092 20092];
beta=1/6;          %定义系数beta
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=1/(beta*detat^2);  
a1=0.5/(beta*detat);
a2=1/(beta*detat);
a3=1/(2*beta)-1;
a4=0.5/beta-1;
a5=detat/2*(0.5/beta-2);
a6=detat*(1-0.5);
a7=0.5*detat;
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
                                        %计算j时刻位移、加速度和速度向量
   y(:,j)=inv(k2)*(p(:,j)+m*(a0*y(:,j-1)+a2*y1(:,j-1)+a3*y2(:,j-1)));   
   y2(:,j)=a0*(y(:,j)-y(:,j-1))-a2*y1(:,j-1)-a3*y2(:,j-1);
   y1(:,j)=y1(:,j-1)+a6*y2(:,j-1)+a7*y2(:,j);
  end
tt2=0:detat:tt;
subplot(2,2,1),plot(tt2,y(1,:))
title('The Newmark Method')
xlabel('t')
ylabel('X1')
subplot(2,2,2),plot(tt2,y(2,:))
title('The Newmark Method')
xlabel('t')
ylabel('X2')
fid=fopen('shuju2.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);