中心差分法求解两自由度体系的振动反应

t0=0;u0=[0 0];v0=[0 0];
all_time=input('输入周期T(1x1)：');
c=input('输入阻尼C(2x2)：');
dt=all_time/10;
m=[2 0;0 1];
k=[6 -2;-2 4];
P0=[0;10];

all_times=all_time+3*dt;
t=[0:dt:all_times;0:dt:all_times];
[mm,nn]=size(t);
u=zeros(size(t));
v=zeros(size(t));
ac=zeros(size(t));
u(:,2)=u0;
v(:,2)=v0;
ac(:,2)=m\(P0-c*v(:,2)-k*u(:,2));
c0=1/(dt^2); %计算积分常数
c1=1/(2*dt);
c2=2*c0;
c3=1/(c2);
u(:,1)=u(:,2)-dt*v(:,2)+((dt^2)*ac(:,2))/2; %求出初始条件u-1项

ek=m*c0+c*c1; %计算等效刚度
a=k-c2*m; %计算方程系数
b=c0*m-c1*c;
P(:,2)=P0; %给出初始荷载条件
eP(:,2)=P(:,2)-a*u(:,2)-b*u(:,1); %计算初始的有效荷载
u(:,3)=ek\eP(:,2); %计算位移u1=u(:,3)
for i=3:nn