碰摩转子系统分岔图

wh--x123456

有哪位大神帮忙看看程序哪错了？为什么不出结果呢？谢谢

1. 主程序：
   
2. clear
   
3. m1=4.5e5;
   
4. m2=3.0e5;
   
5. c=2.1e6;
   
6. k1=1.0e7;
   
7. k2=1.0e7;
   
8. kr=2.5e8;
   
9. delta=0.0025;
   
10. d=0.000:0.0005:0.005;%d是作分叉参数
    
11. e=0.002;
    
12. l=5.0;
    
13. f=0.01;
    
14. alpha=0.0001;
    
15. omega=2*pi*5;
    
16. w0=sqrt(k1/m1);
    
17. w=omega/w0;
    
18. 
    
19. 
    
20. y0=[0.0002;0.0002;0.0002;0.0002];
    
21. for j=1:length(d);
    
22. tt=2*pi/w; %定义步长的
    
23. options=odeset('RelTol',10^-5,'AbsTol',10^-5);
    
24. [t,y]=ode45('fangcheng1',[0:tt/500:100*tt],y0,options,m1,m2,c,k1,k2,kr,delta,d(j),e,l,f,alpha,omega);
    
25. i=20000:100:30000;
    
26. plot(d(j),y(i,1)'.');
    
27. xlabel('d');ylabel('y1');
    
28. hold on;
    
29. end

复制代码

1. 方程程序：
   
2. %四阶龙哥库塔法的程序
   
3. function uu=fangcheng(t,u,flag,m1,m2,c,k1,k2,kr,delta,d,e,l,f,alpha,omega);
   
4. 
   
5. 
   
6. %phi0=0.6/180*pi;         
   
7. %theta0=0.4/180*pi;
   
8. %%%%无量纲化
   
9. w0=sqrt(k1/m1);
   
10. c11=c/m1/w0;
    
11. k11=k1/m1/w0^2;
    
12. k22=k2/m1/w0^2;
    
13. l11=l/delta;
    
14. d11=d/delta;
    
15. e11=e/delta;
    
16. eta=m2/m1;
    
17. w=omega/w0;
    
18. 
    
19. %%%参数无纲量化%%%
    
20. M=1+eta*cos(alpha)^2;
    
21. K=k11+k22*cos(alpha)^2;
    
22. 
    
23. %%合成参数
    
24. r=d11-l11*sin(alpha);
    
25. phi=w*t;         
    
26. theta=w*t;
    
27. 
    
28. A1=e11*w^2*cos(phi);
    
29. A2=eta*r*w^2*cos(theta)*cos(alpha)^2;
    
30. A3=k22*r*cos(theta)*cos(alpha)^2;
    
31. B1=e11*w^2*sin(phi);
    
32. B2=eta*r*w^2*sin(theta)*cos(alpha)^2;
    
33. B3=k22*r*sin(theta)*cos(alpha)^2;
    
34. 
    
35. %%碰摩力
    
36. R=sqrt(u(1)^2+u(3)^2);
    
37. if R>=delta
    
38.     Fpengmo=-(R-1)/R*kr/k1*[1 -f;f 1]*[u(1);u(3)];
    
39. else
    
40.     Fpengmo=[0;0];
    
41. end
    
42. 
    
43. Fx=Fpengmo(1)/M;
    
44. Fy=Fpengmo(2)/M;
    
45. 
    
46. %%%% 进行参数简化后的方程
    
47. uu=zeros(4,1);
    
48. uu(1)=u(2);
    
49. uu(2)=(A1+A2-A3)/M+Fx-c11/M*u(2)-K/M*u(1);
    
50. uu(3)=u(4);
    
51. uu(4)=(B1+B2-B3)/M+Fy-c11/M*u(4)-K/M*u(3);

复制代码

wh--x123456

求做过转子碰摩的大神帮帮忙啊