newmark法 matlab,newmark法和wilson法求解单自由度体系加速度反应谱

fid=fopen('77(0.01-4172).txt','r'); %以下为newmark法求解单自由度体系反应谱

[Acceleration,count]=fscanf(fid,'%g');

dt=0.01;    %时间间隔

d=0.05;     %阻尼比

a=1/2;      %参数1

b=1/4;      %参数2

for i=1:600

t(i)=0.01*i;

w=6.283185/t(i);  %无阻尼自振圆频率

dw=d*w;           %阻尼比乘圆频率

wt=w*dt;          %圆频率乘时间间隔

D=[1+dw*dt+(w*dt)^2*b]^(-1);  %此处应该存在问题

A=[D*(1+2*a*dw*dt-(1-2*b)/2*wt^2-(a-2*b)*d*wt^3),D*dt*(1-(1-2*a)*d*wt-2*(a-2*b)*d^2*wt^2);D/dt*(-wt^2+(a-2*b)/2*wt^4),D*(1-(1-a)*2*dw*dt-(a-b)*wt^2+(a-2*b)*d*wt^3)];

B=[D*dt^2*(-(1/2-b)-(a-2*b)*dw*dt),-1*b*D*dt^2;D*dt*(-(1-a)+(a-2*b)/2*wt^2),-1*a*D*dt];

C=2*dw*abs(Acceleration(1))*dt;      %体系初始时刻绝对加速度

X=[0;Acceleration(1)*dt];            %体系初始时刻位移和速度

for n=2:length(Acceleration)