[求助]NewMark法计算结构响应的问题
我按照书上的NEWMARK原理编了一段MATLAB程序来计算一个3自由度结构的响应，可是感觉有问题，下面是代码，请高手指点！%用newmark法计算一个3自由度结构的响应：
%参数初始化：
clear
clc
N=3;
%N为自由度数。
M=;M=diag(M);
C=;
K=9*;K=diag(K);
d=5;
%间隔时间为1/d秒。
totaltime=6;tt=totaltime*d+2;
%totaltime为外力总持续时间，tt为划分的时间点数，划分精细程度由d决定。
u=zeros(N,tt);
du=zeros(N,tt);
ddu=zeros(N,tt);
%u、du、ddu分别表示位移、速度和加速度。
gama=0.5;
beta=0.25*(0.5+gama^2);
%gama和beta就是γ和β。
detat=1/d;
%按公式计算积分常数bi：
b1=1/(beta*detat^2);
b2=1/(beta*detat);
b3=-1/(2*beta)+1;
b4=gama*b1*detat;
b5=1+gama*b2*detat;
b6=(1+gama*b3-gama)*detat;
EK=K+b1*M+b4*C;
%EK表示有效刚度矩阵。
f=zeros(N,tt);
Ef=zeros(N,tt);
t1=1;
%用t1是为了保证矩阵索引号为整数。
for t=0:1/d:totaltime
%因为detat=0.1，所以t=0:0.1:totalt
f(:,t1)=;
t1=t1+1;
end
t1=1;
for t=0:1/d:totaltime
%计算t+detat时刻的有效载荷：
Ef(:,t1+1)=f(:,t1+1)+M*(b1*u(:,t1)-b2*du(:,t1)-b3*ddu(:,t1))+C*(b4*u(:,t1)-b5*du(:,t1)-b6*ddu(:,t1));
%计算t+detat时刻的位移：
u(:,t1+1)=inv(EK)*Ef(:,t1+1);
%计算t+detat时刻的速度和加速度：
du(:,t1+1)=b4*(u(:,t1+1)-u(:,t1))+b5*du(:,t1)+b6*ddu(:,t1);
ddu(:,t1+1)=b1*(u(:,t1+1)-u(:,t1))+b2*du(:,t1)+b3*ddu(:,t1);
t1=t1+1;
end
t=1:tt;
plot(t,Ef(1,1:tt))
%为什么Ef的曲线不是正弦形式？（位移u的也不是）
grid on
clear
[ 本帖最后由 realyyy 于 2006-8-13 19:42 编辑 ] 这种问题最难解决了，只能自己慢慢查 我的也是，我想是不是激励力激起了一部分自由振动，从而两个波叠加在了一起。我用振型叠加法算出来的是正弦 建议楼主找找动力学书上的经典例子验证一下程序。 有没有非线性动力方程的newmark法 计算程序啊 Newmark方法在数值计算的书中都有现成的代码，可以参考！用代码计算下duffing方法来验证算法！但该算法中的几个参数的确定是有讲究的，希望注意下！
页:
[1]