Newmark法求解两自由度体系的振动反应

clear all;
clc;
t0=0;u0=[0 0];v0=[0 0];
all=input('按格式和顺序，输入周期T、阻尼C、时间步长dt=T/10或T*10、参数α、参数d：如[2.8 0 0.28 0.25 0.5]=')
all_time=all(1);
c=all(2);
dt=all(3); %时间步长
a=all(4);
d=all(5);

m=[2 0;0 1];
k=[6 -2;-2 4];
P0=[0;10];

all_times=all_time+80*dt; %这项只是为了多得几项dt的：t、u、v、ac
t=[0:dt:all_times;0:dt:all_times]; %将时间分步，采用等时间步长；
[mm,nn]=size(t); %计算t的向量长度，得出步数；
u=zeros(size(t)); %设定储存u的矩阵；
v=zeros(size(t)); %设定储存v的矩阵；
ac=zeros(size(t)); %设定储存ac的矩阵；
u(:,1)=u0;
v(:,1)=v0; %赋值向量第1项为u0、v0
if d<0.5 & a<0.25*(0.5+d)^2 %输入参数α和d并判断
disp('不满足稳定条件：d>=0.5,α>=0.25*(0.5+d)^2.请重新输入符合条件的参数')
return
elseif d>=0.5 & a>=0.25*(0.5+d)^2
disp('α和d参数满足条件')
end
ac(:,1)=m\(P0-c*v(:,1)-k*u(:,1)); %求出初始加速度ac0
c0=1/(a*dt^2); %计算积分常数