下面的分岔图是否正确

咕噜噜

前提是他用的最大值法的话，或者根本就是个人编写的程序

octopussheng

呵呵，这倒是，我这里也只是猜测，还是请楼主来为我们答疑吧！

chuandong418

呵呵，循环语句的那个回头我把那个程序找到后给大家传上来，那问下如果以两个参数做为分叉参数的话，改程序如何编写，谢谢了!

octopussheng

你说的这个目前应该还算是世界级的难题吧，呵呵！不知道如何解决！

chuandong418

呵呵，学长误解偶的意思了，我页表达有问题，呵呵就是说如果还有其他的参数的话是不是也是类似的方法来编写？

chuandong418

这种分叉图的思路正确吗？现在有点糊涂，还望指点，非常感谢！

chuandong418

还有你所说的运用最大值法来编写分叉程序，如果用最大值法如何编写？还望指点，谢谢

octopussheng

这种最大值法你可以参考下面的代码，论坛里面其实有一些的，为方便就帖一下吧

主程序：
clear all
global r
t0=[0 150];%积分时间
y0=[1,0,0];
%bifurcation
for r=1:0.05:30 %r的变化精度
[t,y]=ode45('Lorenz',t0,y0);
[Xmax]=getmax(y);
plot(r,Xmax,'b','markersize',1)
hold on
clear Xmax
end
xlabel('r')
ylabel('Xmax')

子程序1：
function dy = Lorenz(t,y)
global r
dy=zeros(3,1);
dy(1)=-10*(y(1)-y(2));
dy(2)=-y(1)*y(3)+r*y(1)-y(2);
dy(3)=y(1)*y(2)-8*y(3)/3;


子程序2
function [Xmax] = getmax(y)
a=length(y);
j=1;
for i=(a-1)/2:a
b=(y(i,1)-y(i-2,1))/2;
c=(y(i,1)+y(i-2,1))/2-y(i-1,1);
if y(i-2,1)<=y(i-1,1)&y(i-1,1)>=y(i,1)&c==0
Xmax(j)=y(i-1,1);
j=j+1;
elseif y(i-2,1)<=y(i-1,1)&y(i-1,1)>=y(i,1)
Xmax(j)=y(i-1,1)-b^2/(4*c);
j=j+1;
end
end

咕噜噜

最大值法论坛里面有编写好的分叉程序，你找一下
两个参数的话只能一个一个参数的来

[ 本帖最后由 咕噜噜 于 2007-8-14 10:28 编辑 ]

chuandong418

现在分岔图形已经画出来了，但是感觉还是不象啊，望大家帮忙!谢谢了1
尤其图形上的点都是花花绿绿的，能不能不那么粗遭哦，感觉周期运动和混沌运动都分不开，望帮忙，谢谢!

[ 本帖最后由 chuandong418 于 2007-8-14 13:17 编辑 ]

无水1324

for  n=1:length(w);
T=2*pi/w(n);
ts=[0:T/100:200*T];
x0=[0.001 0 0.001 0];
[t,X]=ode45('www',ts,x0,options,w(n));
figure(1)
plot(w(n),X(10000:100:end,1),'k.');
hold on
end
xlabel('\fontsize{18}\omega');
ylabel('\fontsize{18}x');
grid

chuandong418

无水学长，我就是按你这个程序去做的，出来的结果就是上面这个图,感觉不象!

chuandong418

另外我使用 b 作为分岔参数,为什么会出现下面的错误,请帮忙@,还需要做些什么修改?
其中参数b取的范围为:2.0:0.1:16.0,出现了如下错误:

??? Input argument "b" is undefined.

Error in ==> www at 17
dx(3)=-2*E*x(3)-x(1)-a*(x(1)^2+x(2)^2)*x(1)-b*(1-1/e)*(x(1)-f*x(2))+u*w^2*cos(w*t);

Error in ==> funfun\private\odearguments at 110
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

octopussheng

把你定义微分方程的程序贴一下吧，估计是你这个程序出问题了！

chuandong418

如下是运行参数b的程序:
function ff=www(t,x,flag,w,b)
ff=zeros(4,1);
g=9800;
E=0.12;
f=0.12;
a=0.5;
%b=10.8;
q=0.16;
f0=25;
u=0.136/q;
w=2.0;
e=sqrt(x(1)^2+x(2)^2);
G=g/((2*pi*f0)^2*q);
dx(1)=x(3);
dx(2)=x(4);
dx(3)=-2*E*x(3)-x(1)-a*(x(1)^2+x(2)^2)*x(1)-b*(1-1/e)*(x(1)-f*x(2))+u*w^2*cos(w*t);
dx(4)=-2*E*x(4)-x(2)-a*(x(1)^2+x(2)^2)*x(2)-b*(1-1/e)*(f*x(1)+x(2))+u*w^2*sin(w*t)-G;
ff=[dx(1);dx(2);dx(3);dx(4)]
===========

b=2.0:0.1:15.0;
options = odeset('RelTol',1e-6,'AbsTol',[1e-6 1e-6 1e-6 1e-6]);
for  n=1:length(b);
T=2*pi/w;
ts=[0:T/100:100*T];
x0=[0.001 0 0.001 0];
[t,X]=ode45('www',ts,x0,options,b(n));
figure(1)
plot(b(n),X(5000:100:10000,1),'.');
xlabel('\fontsize{18}\omega');
ylabel('\fontsize{18}x');grid
hold on
figure(2)
plot(b(n)*30/pi,X(5000:100:10000,2),'.');
xlabel('\fontsize{18}\omega');
ylabel('\fontsize{18}y');grid
hold on
end
运行完后就出现上面发的那些错误，不知道为什么?