请教一个纽马克程序问题

dongwuu

下面是用纽马克法求振动响应的matlab代码, 初始条件u,v,a都为零时，结果良好
可是当v=[4;4;4;4;4;4;4];时，得出的u值特别大，与初始条件u,v,a都为零时相差甚远，这是什么原因？
请教论坛里的朋友了！！！

clc;clear;
j1=76318.474; j2=5255.968; j3=1063.431;
j4=7858.041; j5=821.753; j6=2290.285; j7=11342.453;
k1=6.464e8; k2=6.486e8; k3=5.929e8;
k4=1.981e8; k5=1.980e8; k6=11.595e8; M0=2000;  
J=[j1 0 0 0 0 0 0; 0 j2 0 0 0 0 0; 0 0 j3 0 0 0 0; 0 0 0 j4 0 0 0; ...
   0 0 0 0 j5 0 0; 0 0 0 0 0 j6 0; 0 0 0 0 0 0 j7]; %% J=diag([j1,j2,j3,j4,j5,j6,j7]); %%by ChaChing
K=[k1 -k1 0 0 0 0 0; -k1 k1+k2 -k2 0 0 0 0; 0 -k2 k2+k3 -k3 0 0 0
   0 0 -k3 k3+k4 -k4 0 0; 0 0 0 -k4 k4+k5 -k5 0; 0 0 0 0 -k5 k5+k6 -k6; 0 0 0 0 0 -k6 k6];
degree=length(J); C=zeros(degree,degree);
u(:,1)=zeros(degree,1); v(:,1)=[0;0;0;0;0;0;0]; a(:,1)=zeros(degree,1);
n=10;%时间段数
tend=0.05; dt=tend/n; gaJa=0.5; beta=0.25;
alpha0=1/beta/dt^2; alpha1=gaJa/beta/dt; alpha2=1/beta/dt; alpha3=1/2/beta-1;
alpha4=gaJa/beta-1; alpha5=dt/2*(gaJa/beta-2); alpha6=dt*(1-gaJa); alpha7=gaJa*dt;
Kn=K+alpha0*J+alpha1*C; %%%%形成等效刚度阵
%%%%不同时刻等效荷载Pn,位移u,速度v,加速度a的迭代计算
for i=1:1:n
    P(:,i+1)=[(M0/tend)*dt*i;0;0;0;0;0;(-M0/tend)*dt*i];
    Pn(:,i+1)=P(:,i+1)+J*(alpha0*u(:,i)+alpha2*v(:,i)+alpha3*a(:,i)) ...
        +C*(alpha1*u(:,i)+alpha4*v(:,i)+alpha5*a(:,i));
    u(:,i+1)=Kn\Pn(:,i+1);
    a(:,i+1)=alpha0*(u(:,i+1)-u(:,i))-alpha2*v(:,i)-alpha3*a(:,i);
    v(:,i+1)=v(:,i)+alpha6*a(:,i)+alpha7*a(:,i+1);
end
N=400;   %总时间段数
for j=n+1:1:N  
    P2=[M0;0;0;0;0;0;-M0];
    Pn(:,j+1)=P2+J*(alpha0*u(:,j)+alpha2*v(:,j)+alpha3*a(:,j)) ...
        +C*(alpha1*u(:,j)+alpha4*v(:,j)+alpha5*a(:,j));
    u(:,j+1)=Kn\Pn(:,j+1);
    a(:,j+1)=alpha0*(u(:,j+1)-u(:,j))-alpha2*v(:,j)-alpha3*a(:,j);
    v(:,j+1)=v(:,j)+alpha6*a(:,j)+alpha7*a(:,j+1);
end
%%%%输出
u1=u'

[ 本帖最后由 ChaChing 于 2009-4-9 13:19 编辑 ]

ChaChing

LZ这个不是编程问题了!?
应属专业问题, 待专家并有兴趣者路过