急求：用matlab来画出Lorenz的关于参数r的分叉图程序

lirong75

急求：用matlab来画出Lorenz的关于参数r的分叉图程序老板要求做出来

lirong75

谢谢各位大虾

jumpwolf

1. clear all
   
2. global r
   
3. t0=[0 150];%积分时间
   
4. y0=[1,0,0];
   
5. 
   
6. %bifurcation
   
7. for r=20:0.05:30 %r的变化精度
   
8. [t,y]=ode45('Lorenz',t0,y0);
   
9. [Xmax]=getmax(y);
   
10. 
    
11. plot(r,Xmax,'b','markersize',1)
    
12. hold on
    
13. clear Xmax
    
14. end
    
15. xlabel('r')
    
16. ylabel('Xmax')
    
17. 
    
18. 
    
19. function [Xmax] = getmax(y)
    
20. a=length(y);
    
21. j=1;
    
22. for i=(a-1)/2:a
    
23. 
    
24. b=(y(i,1)-y(i-2,1))/2;
    
25. c=(y(i,1)+y(i-2,1))/2-y(i-1,1);
    
26. 
    
27. if y(i-2,1)<=y(i-1,1)&y(i-1,1)>=y(i,1)&c==0
    
28. Xmax(j)=y(i-1,1);
    
29. j=j+1;
    
30. elseif y(i-2,1)<=y(i-1,1)&y(i-1,1)>=y(i,1)
    
31. Xmax(j)=y(i-1,1)-b^2/(4*c);
    
32. j=j+1;
    
33. end
    
34. end
    
35. 
    
36. function dy = Lorenz(t,y)
    
37. global r
    
38. dy=zeros(3,1);
    
39. dy(1)=-10*(y(1)-y(2));
    
40. dy(2)=-y(1)*y(3)+r*y(1)-y(2);
    
41. dy(3)=y(1)*y(2)-8*y(3)/3;

复制代码

仅供参考！要想图漂亮一点，可以把积分时间和r的变化精度增大，但是计算时间会变长。

lirong75

jumpwolf谢谢你

lmoehg

两点建议
  1 该程序没有考虑过渡过程的剔除，建议每次积分都忽略前面一段足够长时间的解。因为如果是周期解或不动点，这样会得到错误的分岔图。
  2 该程序是用某个分量的局部极大值描绘分岔图，其实没有必要使用两两比较的办法，这样如果数据量大会运行很长的时间。以x分量为例，x的局部极大值就是x的导数为零的点，只需考虑导数值（就是第一个方程的右边）的符号即可，另外应该考虑插值以得到足够的精度。

htwei

楼上的说的很好