Newmark数值分析法matlab程序

clc
clear
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;%c=2m0.05*sqrt(k/m) ;
f0= 100;
t1=5; %力的作用时间
nt=2000; %分2000步完成
dt=0.01; %时间步长
alfa=0.25; %γ=0.25
beta=0.5; %β=0.5

a0=1/alfa/dt/dt; %
a1=beta/alfa/dt; %
a2=1/alfa/dt; %
a3=1/2/alfa-1; %
a4=beta/alfa-1; %
a5=dt/2*(beta/alfa-2); %
a6=dt*(1-beta); %
a7=dtbeta; %
d=zeros(4,nt); %初位移为0
v=zeros(4,nt); % 初速度为0
a=zeros(4,nt); % 初加速度为0
for i=2:nt
t=(i-1)dt;
if (t<t1),
f=[f0sin(4pit/t1);0;0;0]; %f=[300sin(6pit)-50cos(3pit);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));
% t+dt时刻的有效荷载
d(:,i)=ke\fe;%d(:,i)=inv(ke)fe; %求解t+dt时刻的位移
a(:,i)=a0(d(:,i)-d(:,i-1))-a2v(:,i-1)-a3a(:,i-1); %计算t+dt时刻的加速度
v(:,i)=v(:,i-1)+a6a(:,i-1)+a7*a(:,i); %计算t+dt时刻的速度
end

T=(0:dt:19.99); %离散系统dt为采样周期 19.99为终端时间
close all

figure %控制窗口数量
plot(T,d(1,:)) %绘制位移函数图像
title(‘各质点位移总图’) %添加标题为各质点位移总图
xlabel(‘s’) %对x轴进行标注为时间（s）
ylabel(‘m’) %对y轴进行标注为位移（m）
grid %显示画图中的个网线

figure
plot(T,a) %绘制加速度函数图像
title(‘加速度总图’) %添加标题为各质点加速度总图
legend(‘a1’,‘a2’,‘a3’,‘a4’)
xlabel(‘s’) %对x轴进行标注为时间（s）
ylabel(‘m/s^2’) %对y轴进行标注为加速度（m/s2）
grid

figure
plot(T,v) %绘制速度函数图像
title(‘速度总图’) %添加标题为各质点速度总图
legend(‘v1’,‘v2’,‘v3’,‘v4’)
xlabel(‘s’) %对x轴进行标注为时间（s）
ylabel(‘m/s)’) %对y轴进行标注为加速度（m/s）
grid