请教newmark求解动力学方程遇到的问题
下午看了看书，在网上找了找代码，试着用newmark法求解动力学方程，但是结果不对。好像以前在网上有人遇到过相似的问题，请教达人指点！
%用newmark法计算一个3自由度结构的响应：
%参数初始化：
N=3;%N为自由度数。
M=;M=diag(M);
K=;
C=0.*(M+K);
d=10;%间隔时间为1/d秒。
totaltime=20;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就是γ和β。
dt=1/d;
%按公式计算积分常数bi：
alpha0=1/beta/dt^2;
alpha1=gama/beta/dt;
alpha2=1/beta/dt;
alpha3=1/2/beta-1;
alpha4=gama/beta-1;
alpha5=dt/2*(gama/beta-2);
alpha6=dt*(1-gama);
alpha7=gama*dt;
EK=K+alpha0*M+alpha1*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*(alpha0*u(:,t1)+alpha2*du(:,t1)+alpha3*ddu(:,t1))+C*(alpha1*u(:,t1)+alpha4*du(:,t1)+alpha5*ddu(:,t1));
%计算t+detat时刻的位移：
u(:,t1+1)=inv(EK)*Ef(:,t1+1);
%计算t+detat时刻的速度和加速度：
ddu(:,t1+1)=alpha0*(u(:,t1+1)-u(:,t1))-alpha2*du(:,t1)-alpha3*ddu(:,t1);
du(:,t1+1)=u(:,t1)+alpha6*du(:,t1)+alpha7*ddu(:,t1+1);
t1=t1+1;
end
t=0:0.1:20+0.1;
% plot(t,Ef(1,1:tt))；
plot(t,u(1,:))；
grid on
clear
[ 本帖最后由 无水1324 于 2007-7-20 08:41 编辑 ] C=0.*(M+K);
这个阻尼有什么用？
回复 #2 无水1324 的帖子
以前看到交大有个老师的讲义，例子是c＝0，给出了解析解和Rk方法的结果。所以我这里也让c＝0，验证我的代码是否正确 对于每一个时刻你这里面应该要做迭代阿 弱弱的问一句，对于每步计算出来的结果需要进行收敛判断吗？ 本帖最后由 VibInfo 于 2016-5-17 13:22 编辑原帖由 pp786195 于 2009-11-9 14:49 发表
弱弱的问一句，对于每步计算出来的结果需要进行收敛判断吗？
只要满足newmark方法时间步长大小的要求，每步计算出来的结果不需要进行收敛判断
回复 楼主 tianlei2005 的帖子
同意无水给一个合适的C值和合适的时间不长（要满足NEWMARK方法的时间步长要求）应该收敛的
页:
[1]