看图说话——布鲁塞尔问题的分叉图

cam_1980

能不能帮我把不正确的地方指出来？

cam_1980

能不能帮我把不正确的地方指出来？并说明一下理由。我实在不知道错在哪

cam_1980

终于知道问题之所在，我把你给我的程序修改之后运行，图像果然没有毛刺，但是，如果将后面两句
bifdata=YY1(:,end-61:end);————》改为bifdata=YY1；
plot(range,YY1,'k.','markersize',1);
绘制结果与我得出的图像一样毛刺很多。也就是说你开始虽然去除了60个周期的瞬态，在plot之前又舍弃了一部分，只绘制了少力量数据。不知道我分析的对不对。
bifdata=YY1(:,end-61:end);时的分叉图：

bifdata=YY1时分叉图：


[ 本帖最后由 cam_1980 于 2008-8-29 20:54 编辑 ]

无水1324

恩有一 点道理的,

liliangbiao

程序不是死的 ，要根据自己的需要临时调整！我觉得编写程序有一个很重要的问题就是 程序的本身的合理性问题，至于你的那个程序，也包括我给你的那个程序(即我空间中的程序)，我觉得有一点问题就是没有将得到的值作为初值，作为下一次RK45的迭代初值！

无水1324

关于你说的这种初值变化的也是做分岔分析的一个策略，在很多文献中提到过

liliangbiao

是的，这一点我认为很重要！我记得许多文献用一些数学理论能解释为什么得到的值必须作为下一步迭代的初值，会更好些(或者说精确)，但是当时这些文献并没有引起我的注意，因此我没有研读，所以我还没有学会这一点。一旦哪天我再次碰到这样的文献，我仔细研读一下，从数学的角度怎么处理的，我会及时的告诉大家！

cam_1980

通过这个问题，我学到很多东西，非常感谢大家的参与。

liliangbiao

这样的话，我们的目的都达到了不少！但是离最终结束的时间还很长，我们要继续努力哦！！

胡晓宇

取出瞬态的影响，是不是就是去除前面一部分周期的数据，有没有其他的方法可以定量确定目前的状态是稳定状态？不知道 liliangbiao和无水大哥有没有方法？

无水1324

liliangbiao这个问题确实值得深入学习哈，我也在用这个但是总是理所当然的，没有去深究其数学上的意义，有什么进展一定要告知哈

无水1324

简单也是最笨的方法就是，比较2周期内的数据是否吻合，来判断。其他的后面的高手继续

胡晓宇

cam_1980:
能不能把liliangbiao哥给你改的程序发给我一份，我现在做分岔图也遇到毛刺的问题，去掉瞬态也不行，我想用liliangbiao的程序试验一下。谢谢你们两位了
我的信箱：hutsinghua@yahoo.com.cn

liliangbiao

个人主页怎么封闭了，我现在上不去，看着大家对这个挺感兴趣的，找到了这个程序之后，我贴上！

liliangbiao

经不住念叨，个人主页竟然能上了，我吧程序贴上！
你还可以搜索我的个人主页！

% Author: Thomas Lee % E-mail: lixf1979@126.com % Corresponding: School of Mathematics, Physics and Software Engineering, Lanzhou Jiaotong University, Lanzhou 730070, China function dx=duffing(t,X) global F wd; r=0.168; x=X(1); y=X(2); psi=X(3); dx=zeros(3,1); dx(1)=y; dx(2)=-r*y+1/2*x*(1-x^2)+F*sin(psi); dx(3)=wd; clear; global F wd; wd=1.0; range=[0.10:0.0001:0.3]; period=2*pi/wd; % k=0; YY2=[]; step=2*pi/100; %步长。 for F=range y0=[0 0 0]; F k=k+1; % discard the first 60 periodic data; %除去前面60个周期的数据，并将最后的结果作为下一次积分的初值 tspan=[0:step:60*period]; [t,Y]=ode45(@duffing,tspan,y0); y0=Y(end,:); j=1; for i=60:200 tspan=[i*period:step:(i+1)*period]; [t,Y]=ode45(@duffing,tspan,y0); YY1(k,j)=Y(end,1); % get the omega data from every period end j=j+1; %取出每一个周期内的第一个解的最后一个值。 y0=Y(end,:); end end bifdata=YY1(:,end-51:end); plot(range,bifdata,'k.','markersize',1);