连续动力系统的分岔图的画法
本帖最后由 liliangbiao 于 2011-6-1 01:11 编辑以Duffing方程为例，我谈谈我做分岔图的想法和体会。这里没有什么标准的Duffing方程，而是很多种之一，方程如下：
d^2x/dt^2+epsi*dx/dt+a*x+b*x^3=fcos(w*t) where, w=0.8, epsi=0.2,a=1,b=1,f=<10,32>当然 f 可以取到相当大，例如f=<-1000,1000>,但是用不了这么大，如果你做出来了会发现这个系统具有对称性，导致在一定吸引域里面它的数值解也对称，这样导致很多对称的部分(这点被研究过，请参考相关的论文)，于是我们就不妨选取一个小的部分来做研究。
1 频闪法做分岔图的思想：
既然周期力f是一个随着一个周期函数连续变化的，那么就可以采用每个采样周期内打点的方法。注意这里面有两层意思：一层意思是周期函数是连续变化的，另外f也在连续变化。虽然我们的计算机不能做连续的工作，于是我们可以充分的利用计算机的循环功能，来使得f的值取得更为稠密些，比如df=(fmax-fmin)/10000。
2 Poincare 映射(surface of poincare section)的选取
既然要在每个采样周期内打点，那么庞加莱截面技术不得不采用，因为这一技术的几何意义就是告诉我们每个周期内的截面上的被连续轨线横截的情况。庞加莱截面是个非常有用但是很难从一句话中得到含义的子空间。因为不但要考虑轨线的方向性还要考虑轨线在庞加莱截面上的穿越情况。许多做过自治方程(比如Lorenz系统)的，都深有体会，明明是周期一，而取得截面上的点却是两个点。当然，有的截面取得根本不对，导致一片点，在缩小比例的情况下，往往自以为是的认为自己是对的，其实这个是自欺欺人的做法，我会在另文中给出自治微分动力系统的选择截面(surface of poincare section)的方法。那么对于非自治的比如当下讨论的Duffing系统而言，就不用考虑轨线的方向性，这是因为我们取得点是一个周期的情况下，每个点穿越surface of poincare section的情况，因为另次穿越了surface of poincare section但是没有到一个周期，这个点在没有考虑的情况下，就会被忽略掉的(因为在连续动力系统下你不用去考虑)。如果你还不理解我的话，那么就需要补充一点知识了，建议大家到维基网站上看，重点看poincare map 与surface of poincare section的区别和联系。
3 瞬态响应
这个在我的另一个帖子中，给出了详细的说明，建议大家去看连续指数谱的编程思想一贴。不再赘述。
4 RK方法
这是一个经典的解决微分方程不存在解析解方法，建议大家看看这个方法的构造与思想，因为我不建议使用Matlab的，所以对Ode45没感觉，尽管这个是非常大牛的人编写出来的，比你编写的任何的RK方法都要精确、普适，但是由于Matlab对循环太敏感，太慢，我一直不用Matlab编写分岔程序。如果你和我一样，选择C语言或者C++的话，那么就需要看懂这个RK方法的构造了，当然，许多的书上或者网页上都能找到相关的程序供你选择或修改。另外，这些原因也是我共享的Matlab程序中，关于连续系统分岔部分不太精确的原因，因为要调试一个Matlab的分岔程序，需要花费太长时间，要等待太长的时间，所以我更愿意选择一个能在3-10min'里面就能调试的C/C++程序。在你的RK程序中，你的Poincare截面要取得非常小，比如一般都是2*PI/(100*w),这是因为RK计算精度的要求，也就是说一个周期内我利用Poincare算了100次，所以你的循环就是
for(i=1;i<=100;i++)
{
RK(Poincare/100);
}
当然，为了除掉瞬态响应，你要先预算100-1000次，然后再取最后这100次的最后一点，就是你的一个周期内的值。
算法如下
for(j=-1000;j<=100;j++)
{
for(i=1;i<=100;i++)
{
RK(Poincare/100);
}
if(j>0)//前一千次是不要的,为了除掉瞬态响应;保留后100个数值; 这样的话，就相当于每次取100次中的RK最后一个数值,对于每个df而言，被保存了最后100个数值(详细理解见我随后附的数据文件).
fprintf(OUT,"%f,%f\n",df,x1);
}
这里的嵌套循环足以让你不选择Matlab;
5 绘图
数据被保存在txt文件中，你就可以选择任何一个可以调用文本文件的数学软件或者其他绘图软件来画图，以Matlab为例，拷贝到current folder下(我用的是2011b,版本有所不同，目录有所不同)，格式如下
load filename.txt
plot(filename(:,1),filename(:,2),'k.','markersize',1)
6 效果
见图像一二三，可以看到每个df下，系统的运动情况，那么你就可以猜想到在poincare map上对应的点的情况，如果你继续放大，发现自己的周期点不是一个、两个、四个等等，而是一片或者一堆，就说明你计算的不精确，截面选择的不对，RK作用时间不长或者步长选择的不对，当然在C/C++下修改程序非常快，因为你不用等待太长的时间。
如果对我的思想和方法有更深层次的见解的话，可以随时跟帖，我会随时跟进，解决大家的烦恼并且提高我个人的理解和改进我的程序。非常欢迎大家的建议和批评,如果你有思想或者编写的flow chart, flow-process diagram, pseudo codes or codes of C/C++, java etc.程序，请一起附上，望请大家一起修改和改进。另外，请见谅，我不会再对各位编写的程序进行修改或者评价，以免引起误解或者不满，敬请理解 (因为我是来这里交流思想的)。
P-1,周期1的运动(在原图上被放大，点的间隙清楚可见) P-2-4,周期2,4的运动(在原图上被放大，点的间隙清楚可见) RP-1-2,周期1,2的运动(在原图上被放大，点的间隙清楚可见) 本帖最后由 liliangbiao 于 2011-6-1 00:45 编辑
上传的分岔数据，需要的话可以下载，先解压，然后拷贝到你的Matlab工作区内，按照我的提示的方法即可调用画图。另外在这个数据文件中，你会发现，每个df我都在保证在除去了瞬态响应后的，取了最后算的100个点，对于周期一来说，最后100个点都是一个数值，对于周期二来说，最后一百个点都是两组相互交替的数值，等等，这样在图像上，每个df下，你看到的不是一个点，而是100个点。所以混沌区域很稠密，看起来更好看些，当然，你可以选择去做更为细致的图像，比如取f在区间连续变化时，你可以算10000个点df，再者，你可以在保证除掉瞬态响应的情况下，每个df取200-500个点等等，这样会导致你的计算速度放慢或者保存的数据文件过大，不利于初步调试程序和画图，当然你的计算机如果够好的话，比如我用的是i7-990x的cpu, 你可以试试。 楼主的意思是：先得到精确的数值解，然后用频闪法或者选择合适的截面，然后把挑出的点画出该参数对应的点图，通过for循环叠加得到分岔图。这也是分岔图的基本做法吧。对于分岔图，的确用c言要快很多，matlab动不动就几十个小时。
但是，频闪法的周期是怎样确定的，对于多频激励的系统，怎样保证每个周期只取一个点？ 另外，楼主贴出的图在混沌区域，没有很清晰的周期窗口，能否把混沌部分放大看看？
还有，发文章时，很多期刊需要矢量图，但我做的矢量图总是没有JPG的图清晰好看，谁有什么好的意见或者经验？ 回复 7 # kezairenjian 的帖子
数据文件太大，无法上传！需要的话我给你发邮箱！
另外对你的上个回复，我只能说你说对了一半。关于频闪法，你找本书看看吧。一句话说不清楚。至于你说的怎么保持每个周期的最后一个数值，因为R-K方法需要带入你选择的截面来循环，比方说你取得时2*pi/(100*w)循环，那就是R-K迭代到你循环的终点(即循环100次)为止，不就是每个周期取得最后一个么？
楼主，我一直都是用MATLAB做系统的分岔图。我用做一维、二维系统的分岔图方法来画一个三维系统的分岔图，结果总是不对。真不知道是怎么回事.我是用频闪法画分岔图的 对于你所做的思想我不做任何评论，期望你能理解我的这个思想，看看能不能应用到你的所做的系统当中去！ Poincare 映射(surface of poincare section)的选取中，你所说的最后的100个点，如果用一个截面表示，就是着100个点到截面的距离的方程，与截面的交点？是不是这个意思？谢谢？
非常感谢，我是学工科的，能否推荐一本较好的非线性教材？ 大牛能不能推荐一本非线性教材？
我实在不知道哪本教材好，国内的教材呢，一般，国外的教材呢，太厚。有三本教材你可能不能错过，我指的是非线性动力学方面的混沌与分支。
1 Edward Ott, chaos in dynamical systems 国内有卖。99元定价，相对于原版来说便宜得很，并且影印质量不错。
2 K. Alligood, T. Sauer and J. A. Yorke, Chaos: An Introduction to Dynamical Systems。
3 YA Kuznetsov. Elements of Applied Bifurcation Theory, Third Edition国内有翻译版本。 楼主，我想要您matlab绘图的程序，能不能发到我邮箱里？
yanmei922@foxmail.com
谢谢啦