微分方程组的求解发散了，怎么办，求助~

zhengshuhe

我有个5个自由度的常系数非奇次方程组，用ODE45函数求解，得到的结果是发散的，怎么办？
刚度矩阵、质量矩阵、阻尼矩阵、激振力矩阵都是对称的。

zhengshuhe

主程序：
tspan=[0 100];
y0=zeros(10,1);
[t,y]=ode45(@zcx_vfun,tspan,y0);
plot(t,y);grid;
hold on;


函数程序：
function dy=zcx_vfun(t,y)
%系统参数
m_b=650;m_d=3000;m_fd=3000;m_bd=3000;m_f=2100;jd=45/360*6.28;w=200;J_f=4.0*10^3;
F_0=84*10^3;k_bfv=1.2*10^5;c_bfv=1.2*10^2;
k_bfh=1.2*10^5;c_bfh=1.2*10^2;
k_fdv=1.6*10^6;k_bdv=1.6*10^6;
c_fdv=0.87*10^3;c_bdv=0.87*10^3;
k_fdh=1.6*10^6;k_bdh=1.6*10^6;
c_fdh=0.87*10^3;c_bdh=0.87*10^3;
k_fsv=3.2*10^6;k_bsv=3.2*10^6;
c_fsv=11*10^3;c_bsv=11*10^3;
k_fsh=3.2*10^6;k_bsh=3.2*10^6;
c_fsh=11*10^3;c_bsh=11*10^3;
l_1=1.8;l_2=0.8;
%质量矩阵
M=[m_fd,0,0,0,0;0,m_fd,0,0,0;0,0,m_f,0,0;0,0,0,m_b,0;0,0,0,0,J_f];
%刚度矩阵
K=[k_fsv+k_fdv,0,-k_fdv,0,-k_fdv*l_1;0,k_bsv+k_bdv,-k_bdv,0,k_bdv*l_1;-k_fdv,-k_bdv,k_fdv+k_bdv+k_bfv,-k_bfv,k_fdv*l_1-k_bdv*l_1-k_bfv*l_2;0,0,-k_bfv,k_bfv,k_bfv*l_2;-k_fdv*l_1,k_bdv*l_1,k_fdv*l_1-k_bdv*l_1-k_bfv*l_2,k_bfv*l_2,k_fdv*l_1^2+k_bdv*l_1^2+k_bfv*l_2^2];
%阻尼矩阵
C=[c_fsv+c_fdv,0,-c_fdv,0,-c_fdv*l_1;0,c_bsv+c_bdv,-c_bdv,0,c_bdv*l_1;-c_fdv,-c_bdv,c_fdv+c_bdv+c_bfv,-c_bfv,-c_bdv*l_1+c_fdv*l_1-c_bfv*l_2;0,0,-c_bfv,c_bfv,c_bfv*l_2;-c_fdv*l_1,c_bdv*l_1,c_fdv*l_1-c_bdv*l_1-c_bfv*l_2,c_bfv*l_2,c_fdv*l_1^2+c_bdv*l_1^2+c_bfv*l_2^2];
%激振力矩阵
F=[-F_0*cos(jd);-F_0*cos(jd);0;0;0];
%函数定义
x=y(1:5);
dx=y(6:10);
dy=[x;M\(F-K*x-C*dx)];

zhengshuhe

如果直接编写RUNGE-KUTTA 程序而不直接调用ODE函数，会不会好点呢？哪位主任、教授、校长指点下，急盼！

zhengshuhe

激振力矩阵漏了参数：F=[-F*cos(jd)*cos(w*t);-F*cos(jd)*cos(w*t);0;0;0];
但是依然Y矩阵中前5列都为0值（本来应该为位移值按正余弦变化的），后5列向量值有两个发散了，郁闷呢。。。。。高手指点下

zhengshuhe

搞定了！

meiyongyuandeze

ChaChing

zhengshuhe 发表于 2011-3-18 09:13

                               
登录/注册后可看大图

搞定了！

LZ忘了与大伙分享了!:@)