电磁轴承—转子系统的非线性动力学张海燕

shipo507

求助 ，张海燕，电磁轴承—转子系统的非线性动力学。北京工业大学，硕士学位论文：电磁轴承—转子系统的非线性动力学 如何像她说的那样，用matlab分析，得出系统的相图和波形图

shipo507

本人是初学者，请多多指教
程序如下：
function dx=bif(t,y);
dy=zeros(4,1);
omega1=1,omega2=1,p1=-1,p2=1.24,p3=0.75;
beta4=5.9367,beta5=-3.08,q1=-1,q2=1,q3=2,q4=-1.8;
omega=2,f=10;
dy(1)=y(2);
dy(2)=-0.28*y(2)-omega1^2*y(1)-q1*cos(omega*t)*y(1)+(p1*y(1)^3+...
      p2*y(1)*y(3)^2+p3*y(1)^2*y(2)+2*beta4*y(1)*y(3)*y(4)+...
      beta4*y(2)*y(3)^2+2*beta5*y(2)*y(3)*y(4)+3*beta5*y(1)*y(2)^2+...
      beta5*y(1)*y(4))+(q3*y(1)^3+q3*y(1)*y(3)^2+q4*y(1)^2*y(2))*...
      cos(omega*t)-f*sin(omega*t);
dy(3)=y(4);
dy(4)=-0.28*y(4)-omega2^2*y(3)-q2*cos(omega*t)*y(3)+(p2*y(3)^3+p2*y(1)^2*y(3)...
     +p3*y(3)^2*y(4)+2*beta4*y(1)*y(2)*y(3)+beta4*y(1)^2*y(4)+2*beta5*y(1)*...
     y(2)*y(4)+beta5*y(2)^2*y(3)+3*beta5*y(3)*y(4)^2)+(q3*y(3)^3+q3*y(1)^2*y(3)+...
     q4*y(3)^2*y(4))*cos(omega*t)+f*cos(omega*t)
  运行
options=odeset('reltol',1e-4,'abstol',[1e-4,1e-4,1e-5]);
[T,Y]=ode45(@bif,[0,12],[0.1349 0.001 -0.5 0.012],options);
plot(T,Y)
没能得到结果。
现在又放假，同学回去的都差不多了，问了几个也没有对这方面有研究，，

paulhuang

肯定没有收敛，哪有这样的相图？
多算些点

shipo507

我也觉得奇怪，程序我是模仿matlab帮助和参考别人的程序写的，但是运行都不对
现在正在看书，希望自己能尽快弄明白:lol

无水1324

1. function dx=bifff(t,y)
   
2. omega1=1;omega2=1;p1=-1;p2=1.24;p3=0.75;
   
3. beta4=5.9367;beta5=-3.08;q1=-1;q2=1;q3=2;q4=-1.8;
   
4. omega=2;f=10;
   
5. 
   
6. dx=[y(2);
   
7. -0.28*y(2)-omega1^2*y(1)-q1*cos(omega*t)*y(1)+(p1*y(1)^3+...
   
8. p2*y(1)*y(3)^2+p3*y(1)^2*y(2)+2*beta4*y(1)*y(3)*y(4)+...
   
9. beta4*y(2)*y(3)^2+2*beta5*y(2)*y(3)*y(4)+3*beta5*y(1)*y(2)^2+...
   
10. beta5*y(1)*y(4))+(q3*y(1)^3+q3*y(1)*y(3)^2+q4*y(1)^2*y(2))*...
    
11. cos(omega*t)-f*sin(omega*t);
    
12. y(4);
    
13. -0.28*y(4)-omega2^2*y(3)-q2*cos(omega*t)*y(3)+(p2*y(3)^3+p2*y(1)^2*y(3)...
    
14. +p3*y(3)^2*y(4)+2*beta4*y(1)*y(2)*y(3)+beta4*y(1)^2*y(4)+2*beta5*y(1)*...
    
15. y(2)*y(4)+beta5*y(2)^2*y(3)+3*beta5*y(3)*y(4)^2)+(q3*y(3)^3+q3*y(1)^2*y(3)+...
    
16. q4*y(3)^2*y(4))*cos(omega*t)+f*cos(omega*t)];

复制代码

1. 
   
2. clear;close all;clc;
   
3. [t,y]=ode45('bifff',[0,10],[0.1349,0.001,-0.5,0.012]);
   
4. plot(t,y)
   

复制代码

shipo507

你好，谢谢你的帮助
可是我按你的也运行不出来，都已经运行了几个小时都没有出结果。
麻烦你再看看。