转子动力学中matlab用ode45求响应和轴心轨迹问题，请求帮助

一路向前

大家好，我最近在学习转子动力学，用matlab中的ode45求响应和轴心轨迹，遇到问题，求出结果为NaN，求高手指点，多谢。
题目为：
[M1]{u1''}+[c1]{u1'}+w*[J1]{u2'} + [K1]*{u1}={Px(t)}
[M1]{u2‘'}+[c2]{u2'}-w*[J1]{u1'}  +[K1]*{u2}={Py(t)}
[M1]等均为6*6的矩阵。求u1中的某一响应。
我编写得程序为：
function dz=zhuanzi2pan(t,z)
M1
J1
K1  %这些都是根据特定结构求的得，过程略去了
w=500 ;  
Q1=zeros(6,1);   
Q1(2)=m22*10^(-6)*w^2;
C=2*M1+0.001*K1;%按照比例阻尼算的，不知道对不对        

dz=zeros(24,1);
dz(1:6)=z(7:12);
dz(7:12)=M1\Q1*cos(w*t)-w*M1\J1*z(19:24)-M1\C*z(7:12)-M1\K1*z(1:6);
dz(13:18)=z(19:24);
dz(19:24)=M1\Q1*sin(w*t)+w*M1\J1*z(7:12)-M1\C*z(19:24)-M1\K1*z(13:18);


%程序调用
clear,clc
%format long
z0=zeros(24,1);%初始条件
z0(2)=0.1;
z0(14)=0.1;
z0(8)=1;
z0(20)=1;
z0(4)=0.1;
z0(16)=0.1;
z0(10)=1;
z0(22)=1;

tspan=[0:0.01:1];
[t ,z] = ode45(@zhuanzi2pan,tspan,z0);
figure(1)
plot(t,z(:,2))

最终我求得的Z是NaN,求高手指点错误在哪里，谢谢了！

竹语随风

兄弟，我也是研究这个方向的，现在同样在学习转子动力学和MATLAB编程，希望和你有学术交流的机会。

a0302922

程序无误，结果不收敛，可以重选初值点

一路向前

竹语随风 发表于 2013-4-16 22:10

                               
登录/注册后可看大图

兄弟，我也是研究这个方向的，现在同样在学习转子动力学和MATLAB编程，希望和你有学术交流的机会。

嗯，加你好友了，呵呵

一路向前

a0302922 发表于 2013-4-20 14:17

                               
登录/注册后可看大图

程序无误，结果不收敛，可以重选初值点

多谢，没调出来，换别的算法解了

ME!

一路向前 发表于 2013-4-22 19:32

                               
登录/注册后可看大图

多谢，没调出来，换别的算法解了

换什么算法啊，分享下，我用有限元法求刚度和质量，振型叠加法求响应

一路向前

ME! 发表于 2013-7-12 16:13

                               
登录/注册后可看大图

换什么算法啊，分享下，我用有限元法求刚度和质量，振型叠加法求响应

最后用谐波平衡法解的，我研究的故障转子的

ME!

竹语随风 发表于 2013-4-16 22:10

                               
登录/注册后可看大图

兄弟，我也是研究这个方向的，现在同样在学习转子动力学和MATLAB编程，希望和你有学术交流的机会。

请问楼主有没有画过轴心轨迹和瞬态响应啊
下面是我编的，我不知道对不对
http://forum.chinavib.com/thread-127572-1-1.html