请教关于分岔程序
有前辈在吗？gghhjj前辈给的分岔程序，中间部分看不懂，好象道理懂，但具体怎么实现还是没看懂．能否把中间部分给解释一下？我是初学者．我直接运行你的程序for p=linspace(6,13,280);
=ode45('Duffing',10,);
=ode45('Duffing',100,Y(end,:)); ％这一命令是干什摸？
for k=2:length(Y)
f=k-1;
y=1;
if Y(k,1)<0
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
else
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
end
if y<0
Z=; ％这里运行有错误，是y=吗？为什末？
end
end
end
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条件下，也就是说虽然判断结果一样可是成立条件却是不一样的 我觉得做分岔图程序的思路主要分一下的几块：
（1）确定分岔参数的取值范围，并给出变化步长
（2）做一个循环，求解微分方程
（3）可以另做一个判断大小的函数如getmax()，从微分方程的解向量中取出一个极值
（4）将该极值点与相应的分岔参数在图中画出来
大致的方法就是这样，希望这些能够帮助楼主理解程序吧
关于分岔程序
谢谢各位前辈！我的方程是duffing方程：
dx1/dt=x2
dx2/dt=-a*x1-b*x1^3-eta*x2+gamma*cos(omega*t)
上面发的是gghhjj前辈给的程序．
能否把我的方程关于omega或者gamma的分岔程序给写一下．我自己初学，看了半天也不懂．万分感谢各位！！:loveliness: 搜索论坛，我们写过很多这样的程序！
回复 #5 无水1324 的帖子
不好意思，我搜索了，好象没几条．关键是中间的实现过程看不懂，又不知道怎么办．如果我的问题程序能通，我再慢慢琢磨，我想这样能快些．．．．:loveliness: ．多次麻烦，不好意思．[ 本帖最后由 mjhzhjg 于 2007-7-12 21:42 编辑 ] function xdot=duffing(t,x,flag,omega)
eta= ; a= ; b= ;gamma= ;
xdot=
clear;clc;close all;
omega=1:0.01:2;
for h=1:length(omega)
T=2*pi/omega(h);
=ode45('duffing',,,[],omega(h));
plot(omega(h),x(20000:100:end,2),'k.');hold on
end 其中的参数你自己修改、给定
回复 #8 无水1324 的帖子
:lol :'( 太感谢了！！我得认真琢磨，通了以后告诉你．但愿一切顺利！
顺便问一下《理论力学的计算机模拟 这本书是哪个出版社出的？还是前辈告诉别人时我搜索到的．再次谢谢！！！:lol 我的邮箱是csyd5053@yahoo.com.cn
有时间给我发邮件，我再告诉你吧
回复 #10 无水1324 的帖子
好．谢谢前辈百忙中指点．祝愿前辈工作、学业更上一层楼！！:victory:回复 #9 mechanic05 的帖子
清华大学出版社彭芳麟 管靖 胡静 卢圣治编
回复 #8 无水1324 的帖子
前辈，我把参数都输进去了，可是不出结果，老死机．为什么？？程序如下：
function xdot=ddd(t,x,flag,omega)
eta= 0.1; a=-48.704 ; b= 24.35*10^4;gamma=1 ;
xdot=;
clear;clc;close all;
omega=1:0.01:2;
for h=1:length(omega)
T=2*pi/omega(h);
=ode45('ddd',,,[],omega(h));
plot(omega(h),x(20000:100:end,2),'k.');hold on
end
:@(
[ 本帖最后由 无水1324 于 2007-7-11 18:14 编辑 ]
回复 #14 mechanic05 的帖子
你确定一下方程是否有问题。程序没有什么问题了。
我修改了一下。
[ 本帖最后由 无水1324 于 2007-7-11 18:16 编辑 ] 我这里可以算出来，不是死机，而是计算的时间比较长。
clear;clc;close all;
omega=1:0.01:2;
for h=1:length(omega)
T=2*pi/omega(h);
=ode45('ddd',,,[],omega(h));
plot(omega(h),x(10000:100:end,2),'k.');hold on
end
%————————————————————————————————