求助：Newmark算法

zhzq626

各位网友：我的课题是做动力学方面的研究，系统的模型已经建出来了，是20X20的线性时变系统，接下来就是对模型的求解。近期看了一些相关的文献，这种问题的求解方法好多前人用Newmark方法，我本人之前没有遇到过类似的问题，请大家给予一些指导意见：这种方程的求解用什么方法较好；什么文献对newmark方法的阐述较为清晰。

期待对此熟悉的各位网友不吝赐教，多谢！

[ 本帖最后由 无水1324 于 2008-5-27 07:54 编辑 ]

zhzq626

怎么大家只是进来逛逛就走呀？请多多指教。

hao1982

NEWMARK资料比较多，结构动力学书上就有介绍，不难

zhzq626

谢谢 hao1982  的回复。

zhzq626

%%  求解线形系统%% dx^2+x=0
%% dx(0)=0  x(0)=1

clear;clc;
format long

M=[1];
C=[0];
K=[1];
P=[0];
beita_Newmark=1/4;
gama_Newmark=1/2;
derta_t=0.1;  
i_Newmark=75;
a0_Newmark=1/(beita_Newmark*(derta_t)^2);
a1_Newmark=gama_Newmark/(beita_Newmark*derta_t);
a2_Newmark=1/(beita_Newmark*derta_t);
a3_Newmark=1/(2*beita_Newmark)-1;
a4_Newmark=gama_Newmark/beita_Newmark-1;
a5_Newmark=derta_t*(gama_Newmark/beita_Newmark-2)/2;
a6_Newmark=derta_t*(1-gama_Newmark);
a7_Newmark=derta_t*gama_Newmark;
Xa_0_Newmark(1,1)=1;  
Xa_1_Newmark(1,1)=0;
Xa_2_Newmark(1,1)=0;  
Xb_0_Newmark(1,1)=0;  
Xb_1_Newmark(1,1)=0;  
Xb_2_Newmark(1,1)=0;   
K_Newmark(1,1)=0;   
P_Newmark(1,1)=0;
X0(1,i_Newmark+1)=0;
X1(1,i_Newmark+1)=0;
X2(1,i_Newmark+1)=0;
TT(1,i_Newmark+1)=0;
Xa_2_Newmark=inv(M)*P
X0=Xa_0_Newmark;  
X1=Xa_1_Newmark;   
X2=Xa_2_Newmark;

for ii=1:i_Newmark
    t=ii*derta_t;
    TT(1,ii+1)=t;
K_Newmark=K+a0_Newmark*M+a1_Newmark*C;   
P_Newmark=P+M*(a0_Newmark*Xa_0_Newmark+a2_Newmark*Xa_1_Newmark+a3_Newmark*Xa_2_Newmark)+C*a1_Newmark*Xa_0_Newmark+a4_Newmark*Xa_1_Newmark+a5_Newmark*Xa_2_Newmark);

Xb_0_Newmark=inv(K_Newmark)*P_Newmark;  
Xb_2_Newmark=a0_Newmark*(Xb_0_Newmark-Xa_0_Newmark)-a2_Newmark*Xa_1_Newmark-a3_Newmark*Xa_2_Newmark;  
Xb_1_Newmark=a7_Newmark*Xb_2_Newmark+Xa_1_Newmark+a6_Newmark*Xa_2_Newmark;  
X0(1,ii+1)=Xb_0_Newmark;
X1(1,ii+1)=Xb_1_Newmark;
X2(1,ii+1)=Xb_2_Newmark;

Xa_0_Newmark=Xb_0_Newmark;
Xa_1_Newmark=Xb_1_Newmark;
Xa_2_Newmark=Xb_2_Newmark;
end

plot(TT(1,:),X0(1,:),'LineWidth',1)

zhzq626

%%  求解非线性系统
%% dx^2+(100+1000*x^2)*x=0
%% dx(0)=60  x(0)=0

[ 本帖最后由 zhzq626 于 2008-6-8 11:24 编辑 ]