请教关于分岔程序

mechanic05

前辈，我的程序总点数有２００００吧，怎摸还是有问题？:@L :@L :@L

无水1324

你的取点应该小于20000才能够画出来图哈

octopussheng

你好，算的时候a就是取的a= -48.704，我把代码贴出来，你看看那些需要改的告诉我：
function xdot=dddd(t,x,flag,omega)
eta=0.1 ; a= -48.704; b= 24.35e4;gamma= 1;
xdot=[x(2);a*x(1)-b*x(1)^3-eta*x(2)+gamma*cos(omega*t)];


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

无水1324

程序，没有问题，可能你的参数不处于混沌区域

mechanic05

谢谢！
这是我画出的分岔图．请指教！其中程序为
function xdot=dddd(t,x,flag,omega)
eta=0.1 ; a= -48.704; b= 24.35e4;gamma= 1;
xdot=[x(2);-a*x(1)-b*x(1)^3-eta*x(2)+gamma*cos(omega*t)];


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

mechanic05

就是xdot=[x(2);-a*x(1)-b*x(1)^3-eta*x(2)+gamma*cos(omega*t)];中的a前面有负号的．对于我画出的图，要去掉瞬态影响该如何改动呢？请前辈指点？

咕噜噜

要去掉瞬态主要还是把前面的的一些不稳定点去掉，理解了意思你应该可以改了吧
plot(omega(h),x(10000:100:end,2),'k.');hold on
这一项是不是，主要还是改变一下10000，增大

[ 本帖最后由 咕噜噜 于 2007-7-16 15:39 编辑 ]

无水1324

plot(omega(h),x(20000:100:end,2),'k.');hold on

mechanic05

前辈好！
我用这个程序
clear;clc;close all;
omega=0.1:0.01:4.1;
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(30000:100:end,2),'k.');hold on
end

算时为什摸老出现这个问题：
??? Error using ==> plot:@)
Vectors must be the same lengths.

Error in ==> wushui2 at 6
plot(omega(h),x(30000:100:end,2),'k.');hold on

mechanic05

前辈好！
我用这个程序
clear;clc;close all;
omega=0.1:0.01:4.1;
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(30000:100:end,2),'k.');hold on
end

算时为什摸老出现这个问题：
??? Error using ==> plot
Vectors must be the same lengths.

Error in ==> wushui2 at 6
plot(omega(h),x(30000:100:end,2),'k.');hold on

无水1324

[0:T/100:120*T]总共才计算12000个点，你还从30000画图，肯定是没有了

mechanic05

:@L :@L :@L
刚才发现弄错了！已经改正了．正在算．
clear;clc;close all;
omega=0.1:0.01:4.1;
for h=1:length(omega)
T=2*pi/omega(h);
[t,x]=ode45('ddd',[0:T/100:300*T],[0 0],[],omega(h));
plot(omega(h),x(30000:100:end,2),'k.');hold on
end
不知道这样改，结果会怎样？！:@o

mechanic05

天那！！:'(
算了一晚上还在算！也太慢!

无水1324

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

咕噜噜

我这里算了很长时间也没出来结果