怎么样才能绘制出最准确、最精确的分岔图？？

yufeiyfyf

        根据论坛里相关的总结帖和分岔图程序（见octopussheng：【总结帖】分岔图绘制不同方法的总结、比较http://forum.vibunion.com/thread-60104-1-1.html），说对系统单参数分岔图的计算主要有这几种方法：最大值法（用getmax函数进行取点）、Poincare截面法（程序里都是用x=y来截取）和闪频发等，根据这几天的仿真和liliangbiao先生的说法，我也感觉这几种方法都不太精确，liliangbiao先生的说：因为选择Poincare截面有讲究，一般不选取最大值或者其他值，因为这些极大值可能导致图像的不连续性，因为Rung-kutta方法本身的缺陷，可能导致在某个区域内无解或者解的不精确性，因此一般选取的截面为一个平衡点处或者某一个特殊的截面，所以我认为你总结的程序不精确！
          我觉得他写的程序大概都是根据关于Rossler的分叉改成的。现在给出Rossler的分叉程序：
%这是Rossler系统的分岔程序！

1. function u =rosser(t,x)
   
2. global  c
   
3. a=0.2; b=0.2;
   
4. x0=[0.5,0.4,0.3]';
   
5. u =[-x(2)-x(3), x(1)+a*x(2), b+x(3)*(x(1)-c)]';

复制代码

1. clear all
   
2. global c
   
3. zhang=[];
   
4. M=[2.5:0.001:6.5];
   
5. counter=1;
   
6. for counter=1:length(M)
   
7.        c=M(counter);
   
8.     i=2;
   
9. xmax=0;
   
10. xmaxold=0;
    
11. frmdata=[];
    
12. error=0;
    
13. tspan=[0 300];
    
14. var=1;
    
15. y0=[0.5;0.3;0.2];
    
16. [t,x]=ode45(@rosser,tspan,y0);
    
17. while i < size(x,1)
    
18.         if x(i-1,var) < x(i,var) & x(i+1,var) <= x(i,var)
    
19.             xmax=x(i,var);
    
20.             if xmaxold ~= 0
    
21.                 frmdata=[frmdata ; xmax xmaxold];
    
22.             end
    
23.                 xmaxold=xmax;
    
24.         end
    
25.         i=i+1;
    
26. end
    
27. r= length(frmdata)-20:length(frmdata);
    
28. fradata1(1,r)=frmdata(r,1);
    
29. zhang=[zhang;fradata1];
    
30. end
    
31. plot(M,zhang,'k.','markersize',1);
    
32.         xlabel(sprintf('c'));
    
33. ylabel(sprintf('x'))

复制代码

仿真后，得出：

，效果很好。


现在问题是，我将这段程序用于Lorenz系统的话，matlab运行就出现问题，现在给出Lorenz的分叉程序：

1. function u =lorenz(t,x)
   
2. global SIGMA;
   
3. global R;
   
4. global BETA;
   
5. SIGMA = 10;
   
6. BETA = 8/3;
   
7. x0=[0,1,0]';
   
8. u =[SIGMA*(x(2)-x(1)), R*x(1)-x(2)-x(1)*x(3), x(1)*x(2)-BETA*x(3)]';

复制代码

1. clear all
   
2. global SIGMA;
   
3. global R;
   
4. global BETA;
   
5. SIGMA = 10;
   
6. BETA = 8/3;
   
7. zhang=[];
   
8. M=[11:0.5:400];
   
9. counter=1;
   
10. for counter=1:length(M)
    
11.        R=M(counter);
    
12.     i=2;
    
13. xmax=0;
    
14. xmaxold=0;
    
15. frmdata=[];
    
16. error=0;
    
17. tspan=[0:0.5:200];
    
18. var=1;
    
19. y0=[0;1;0];
    
20. [t,x]=ode45(@lorenz,tspan,y0);
    
21. while i < size(x,1)
    
22.         if x(i-1,var) < x(i,var) & x(i+1,var) <= x(i,var)
    
23.             xmax=x(i,var);
    
24.             if xmaxold ~= 0
    
25.                 frmdata=[frmdata ; xmax xmaxold];
    
26.             end
    
27.                 xmaxold=xmax;
    
28.         end
    
29.         i=i+1;
    
30. end
    
31. r= length(frmdata)-20:length(frmdata);
    
32. fradata1(1,r)=frmdata(r,1);
    
33. zhang=[zhang:fradata1];
    
34. end
    
35. plot(M,zhang,'k.','markersize',1);
    
36.         xlabel(sprintf('R'));
    
37. ylabel(sprintf('x'))

复制代码

运行后，matlab出现错误，??? Subscript indices must either be real positive integers or logicals.
Error in ==> lorenz_ext_fct at 33
fradata1(1,r)=frmdata(r,1);
看上面程序标有红色的代码，r= length(frmdata)-20:length(frmdata);是将frmdata最后的21个数赋值给r，fradata1(1,r)=frmdata(r,1);这是数组之间赋值。不知道我这样理解可准确。由于出现错误，我修改了20，改为r= length(frmdata)-8:length(frmdata);
运行后出现错误。??? Error using ==> plot
Vectors must be the same lengths.
总之，这段程序耗费我一天时间还没有头绪，不知道谁matlab厉害，能否指点迷津，当然，如果liliangbiao先生在的话，希望能给予指导，谢谢！！欢迎大家讨论。

zhao522

回复 1 # yufeiyfyf 的帖子

这个没研究过，分叉图             知道，但是没深入研究过。

yufeiyfyf

不知道为什么换一个混沌系统，这段程序就不能运行了，看样子还是理论不够扎实啊

meiyongyuandeze

想上传个自己画的分岔图，看图片一直不能上传啊，郁闷，打击热情

yufeiyfyf

回复 4 # meiyongyuandeze 的帖子

你的分岔图是用什么方法画的？讨论下啊。上传图：【高级模式】-图片-上传

meiyongyuandeze

等换个时间，上传图片吧，可能是系统的问题吧。昨天还可以的

yufeiyfyf

在Rossler的分叉程序中将r= length(frmdata)-20:length(frmdata);换成r= length(frmdata)-10:length(frmdata);得出的图一样，看样子与这句代码没有关系。不知道为什么换了个混沌系统，结果就会有错呢？谁能解疑啊？

yufeiyfyf

liliangbiao先生 解释下啊  

雨人

我说说自己的一点看法
上述求取分岔图的程序是基于时域信号的极大值方法，也就是求取仿真时间段内，对应于某一个分岔参数的、某一状态变量的时域信号上的所有极大值。你程序中
r= length(frmdata)-20:length(frmdata);

20是怎么来的!! 由于对应于某一个分岔参数，在一段时间内，时域信号上的极大值可能没有20个，如果是周期运动的话，这取决于你的仿真时间。也就是说，在某些分岔参数下，极大值本来就没有这么多，你那个存储极大值的数组是动态大小的，如果在去最后21个数据的话就出现数组索引超出边界了。
对于周期运动，由于对应的频率成分比较简单，如果选择10个极大值的话，也能反映出简单周期运动的峰值信息，但是如果系统的动力学行为逐渐失稳，慢慢进入混沌区间。此时如果对应每个分岔参数，在选择10个数据来估计系统的动力学行为，个人感觉就有点少了。这也是为什源程序用最后20个数据。
以上是个人意见，有不对的地方还请大家商讨

yufeiyfyf

回复 9 # 雨人 的帖子

那为什么换个系统，该程序就不能运行了呢？比如Lorenz系统

yufeiyfyf

每天顶起，直到满意答案出来。

liliangbiao

报道！
呵呵，关于自治系统画分岔图无非是选择一个合适的Poincare截面,我之所以把这个程序贴上，是因为我当时觉得也是有问题的，但是不知道错在哪里。所以供有兴趣的学者们讨论，但是大家好像都一直依赖别人现成的程序，而不是动手编写(甚至就是套用而已，换个系统，自己能用就成)，我觉得这点很不好！(个人见解而已)。我不知道什么样的分岔图是最准确的、最精确的。

yufeiyfyf

回复 12 # liliangbiao 的帖子

哎，动手编写，matlab学的一般，这怎样才算好？

雨人

回复 12 # liliangbiao 的帖子

liliangbiao先生好
个人感觉，分岔图只是定性的反应系统的动力学行为，目的是对应每一个分岔参数找到合理庞加莱界面，来反映系统的吸引子。由于系统是自治的，庞加莱截面的选取不能使用频闪法。如果对低于三维系统使用函数F(X,Y,Z)=AX+BY+CZ做截面来反映系统吸引子的话，截面对应每一个分岔参数的方向性如果判断呢？
上面的程序使用极值或者使用最大值法，也就是对动力学系统仿真得到的时间历程，从峰值的不同上来反映其频率不同，进一步推断系统是做几周期运动，如果系统的峰值比较多，也就认为时间历程对应的频率成分较多，系统的动力学行为较复杂。这样的话，仅仅从系统时间历程的复杂性上考虑其动力学行为是否失稳，进入混沌区是否严谨。还请liliangbiao先生指导一下。对于直接使用截面法做分岔图，截面的方向如何确定呢？

yufeiyfyf

回复 14 # 雨人 的帖子

同样的疑问啊