请教关于分岔程序

mechanic05

有前辈在吗？gghhjj前辈给的分岔程序，中间部分看不懂，好象道理懂，但具体怎么实现还是没看懂．能否把中间部分给解释一下？我是初学者．我直接运行你的程序

1. for p=linspace(6,13,280);
   
2.        [T,Y]=ode45('Duffing',10,[3;0;p]);
   
3.        [T,Y]=ode45('Duffing',100,Y(end,:));　　％这一命令是干什摸？
   
4.        for k=2:length(Y)
   
5.           f=k-1;
   
6.           y=1;
   
7.           if Y(k,1)<0
   
8.              if Y(f,1)>0
   
9.                y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));　　％这一句和下面判断句结果为什末一样？
   
10.              end
    
11.           else
    
12.              if Y(f,1)<0
    
13.                y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));
    
14.              end
    
15.           end
    
16.           if y<0
    
17.              Z=[Z p+y*i];　％这里运行有错误，是y=[y p+y*i]吗？为什末？
    
18. end
    
19.        end
    
20.    end
    
21. plot(b,y)

复制代码

也不通，请指教．非常感谢！！看了半天也不懂．望不吝指教．万分感谢！！！！:@)

咕噜噜

具体你的题目我不清楚，但是从程序来看  if Y(f,1)>0
               y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));　　
             end
这一部分是在条件Y(k,1)<0下，而第二个重复部分是在条件Y(k,1)>0条件下，也就是说虽然判断结果一样可是成立条件却是不一样的

octopussheng

我觉得做分岔图程序的思路主要分一下的几块：
（1）确定分岔参数的取值范围，并给出变化步长
（2）做一个循环，求解微分方程
（3）可以另做一个判断大小的函数如getmax()，从微分方程的解向量中取出一个极值
（4）将该极值点与相应的分岔参数在图中画出来

大致的方法就是这样，希望这些能够帮助楼主理解程序吧

mechanic05

谢谢各位前辈！
我的方程是duffing方程：
dx1/dt=x2
dx2/dt=-a*x1-b*x1^3-eta*x2+gamma*cos(omega*t)
上面发的是gghhjj前辈给的程序．
能否把我的方程关于omega或者gamma的分岔程序给写一下．我自己初学，看了半天也不懂．万分感谢各位！！:loveliness:

无水1324

搜索论坛，我们写过很多这样的程序！

mechanic05

不好意思，我搜索了，好象没几条．关键是中间的实现过程看不懂，又不知道怎么办．如果我的问题程序能通，我再慢慢琢磨，我想这样能快些．．．．:loveliness: ．多次麻烦，不好意思．

[ 本帖最后由 mjhzhjg 于 2007-7-12 21:42 编辑 ]

无水1324

function xdot=duffing(t,x,flag,omega)
eta= ; a= ; b= ;gamma= ;
xdot=[x(2);a*x(1)-b*x(1)^3-eta*x(2)+gamma*cos(omega*t)]

clear;clc;close all;
omega=1:0.01:2;
for h=1:length(omega)
T=2*pi/omega(h);
[t,x]=ode45('duffing',[0:T/100:220*T],[0 0],[],omega(h));
plot(omega(h),x(20000:100:end,2),'k.');hold on
end

无水1324

其中的参数你自己修改、给定

mechanic05

:lol :'( 太感谢了！！
我得认真琢磨，通了以后告诉你．但愿一切顺利！
顺便问一下《理论力学的计算机模拟　这本书是哪个出版社出的？还是前辈告诉别人时我搜索到的．再次谢谢！！！:lol

无水1324

我的邮箱是csyd5053@yahoo.com.cn
有时间给我发邮件，我再告诉你吧

mechanic05

好．谢谢前辈百忙中指点．祝愿前辈工作、学业更上一层楼！！:victory:

咕噜噜

清华大学出版社
彭芳麟 管靖 胡静 卢圣治编

mechanic05

前辈，我把参数都输进去了，可是不出结果，老死机．为什么？？
程序如下：

function xdot=ddd(t,x,flag,omega)
eta= 0.1; a=-48.704 ; b= 24.35*10^4;gamma=1 ;
xdot=[x(2);-a*x(1)-b*x(1)^3-eta*x(2)+gamma*cos(omega*t)];

clear;clc;close all;
omega=1:0.01:2;
for h=1:length(omega)
T=2*pi/omega(h);
[t,x]=ode45('ddd',[0:T/100:220*T],[0 0],[],omega(h));
plot(omega(h),x(20000:100:end,2),'k.');hold on
end
:@(

[ 本帖最后由 无水1324 于 2007-7-11 18:14 编辑 ]

无水1324

你确定一下方程是否有问题。
程序没有什么问题了。
我修改了一下。

[ 本帖最后由 无水1324 于 2007-7-11 18:16 编辑 ]

无水1324

我这里可以算出来，不是死机，而是计算的时间比较长。

clear;clc;close all;
omega=1:0.01:2;
for h=1:length(omega)
T=2*pi/omega(h);
[t,x]=ode45('ddd',[0:T/100:120*T],[0 0],[],omega(h));
plot(omega(h),x(10000:100:end,2),'k.');hold on
end

%————————————————————————————————