|
|
马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
本帖最后由 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.程序，请一起附上，望请大家一起修改和改进。另外，请见谅，我不会再对各位编写的程序进行修改或者评价，以免引起误解或者不满，敬请理解 (因为我是来这里交流思想的)。
|
-
D1
|
[复制链接]
发表于 2011-6-1 00:12
提升卡
变色卡
显身卡
楼主
