求救，matlab单自由度的时程分析程序
如题，自己按照Newmark-beta法编制了一个单自由度的时程分析程序，但算出来的结果和sap2000的不符，而且反应加速度值很大，不知道问题出在哪里，源程序如下：clear;clc;
% 结构模型初始参数----------------------------------------------------------
m=3e3; %质量(单位：kg)
k=1e6; %刚度（(单位：N/m)）
kesai=0.05; %阻尼比取0.05
c=2*kesai*sqrt(k*m); %阻尼系数
% 读取地震波数据------------------------------------------------------------
acc=textread('D:\研究生课题资料\科研课题资料\处理后的smc文件\51WCW_90_chnua370295.smc_090501.a','%f','headerlines',56);
PGA_Max=max(abs(acc)) %最大地面加速度绝对值
% Newmark-beta法的基本参数--------------------------------------------------
beta=1/6;gama=0.5; %按线性加速度法计算更接近真实结果，故取此组参数
dt=0.02; %地震加速度时程波记录时间间隔
b1=1/(beta*dt^2); b2=1/(beta*dt); b3=1-1/(2*beta); %计算参数
b4=gama/(beta*dt); b5=gama/beta-1; b6=(1-gama/(2*beta))*dt;
ke=k+m*b1+c*b4; %等效刚度
% 设定结构初始状态为零,生成向量空间存储计算值---------------------------------
u=zeros(100/dt,1); v=zeros(100/dt,1); a=zeros(100/dt,1);
% Newmark-beta法的主计算程序------------------------------------------------
for n=2:100/dt
fe=-m*acc(n)+*m+*c;%等效荷载
u(n)=fe/ke;
a(n)=b1*-b2*v(n-1)+b3*a(n-1);
v(n)=b4*-b5*v(n-1)+b6*a(n-1);
end
% 绘制结构在地震作用下的位移、速度、加速度时程曲线-----------------------------
subplot(3,1,1)
t=(0:length(a)-1)*dt;
plot(t,a) %加速度时程曲线
Acc_Max=max(abs(a))
title('Earthquake Response Curve Of Station 51WCW-90','fontsize',15)
ylabel('Acc(cm/s^2)','fontsize',12)
subplot(3,1,2)
plot(t,v) %速度时程曲线
Vel_Max=max(abs(v))
ylabel('Vel(cm/s)','fontsize',12)
subplot(3,1,3)
plot(t,u) %位移时程曲线
Dis_Max=max(abs(u))
xlabel('Time/s','fontsize',12)
ylabel('Dis/cm','fontsize',12)
% End---程序结束-------------
小弟初次发贴，恳请达人们帮分析一下，不胜感激！
其中的循环部分是根据结构动力学书上的写的，感觉问题就出在那部分了，请高人们指点一下 急啊，请达人们帮下忙 附地震波数据文件
原文件.a格式的，传不上来，所以改成txt格式了。期待达人们的解答
页:
[1]