采用Newmark-β法求解振动方程——MATLAB

声明：该博文仅用于matlab学习！
以四层刚架为例进行计算分析。某四层框架结构在顶部受一个简谐荷载的作用：

力的作用时间t1=5.0s，F0=100N。试计算响应的时间为20s（计算分2000步完成）。阻尼矩阵由Rayleigh阻尼构造。

MATLAB程序编写如下：
m=[1,2,3,4];
m=diag(m);
k= [800 -800 0 0; -800 2400 -1600 0; 0 -1600 4800 -3200; 0 0 -3200 8000];
c=0.05m+0.02k;
f0=100;
t1=5;
nt=2000;
dt=0.01;
gama=0.5;
beta=0.25;
a0=1/ beta/dt/dt;
a1=gama/ beta /dt;
a2=1/beta/dt;
a3=1/2/ beta -1;
a4=gama/beta-1;
a5=dt/2*(gama/ beta-2);
a6=dt*(1-gama);
a7=dt*gama;

d=zeros(4,nt+1);
v=zeros(4,nt+1);
a=zeros(4,nt+1);
tt=zeros(1,nt+1);

for i=2:(nt+1)
tt(1,1)=1dt;
tt(1,i)=idt;
t=(i-1)dt;
if (t<t1),
f=[f0sin(4pit/t1);0;0;0];
else
f=[0;0;0;0];
end

ke=k+a0m+a1c;
fe=f+m*(a0d(:,i-1)+a2v(:,i-1)+a3a(:,i-1))+c(a1d(:,i-1)+a4v(:,i-1)+a5a(:,i-1));
d(:,i)=inv(ke)fe;
a(:,i)=a0(d(:,i)-d(:,i-1))-a2v(:,i-1)-a3a(:,i-1);
v(:,i)=v(:,i-1)+a6a(:,i-1)+a7*a(:,i);
end
figure;
plot(tt(1,:),d(1,:))
ylabel(‘位移/[m]’,‘fontsize’,20);
grid on;
xlabel(‘时间/[s]’,‘fontsize’,20);
title(‘顶层楼板位移响应时程’,‘fontsize’,25);