求下面方程的分岔图

siyaoming

x''+a*(1+b*x^2)*x'+w0^2(1+r*x^2)*x=fcos(omega*t)
a b w0 r f omega 是参数；要求画出以f为横坐标的以x或x'为纵坐标的分岔图，也就是f为自变量，其他参数设好定值就可以了。我的程序是这样的：
clc;
w0=1;a=-5;b=-0.1;r=1.1;omega=0.9;T=2*pi;
f=0:0.1:10;tspan=0:T/25:1000;
y1=[];y2=[];G=[];
for i=1:length(f)
    [t,y]=ode45(@federer,tspan,[0.1,-0.2]',[],w0,a,b,r,f(i),omega);n=length(y(:,2));
y1=[y1,y(1000:2000:end,1)'];
y2=[y2,y(1000:100:end,2)'];m=length(y1);
for j=1:m/i   
G=[G,f(i)];   
end;p=length(G);
N=length(f);
Over_to_step=N-i;
end
figure(1)
plot(G,y1,'k.','markersize',3)
xlabel(' ');
ylabel(' ');
title(' ');   


function dx=federer(t,x,w0,a,b,r,f,omega)
dx=[x(2);-w0^2*x(1)-a*x(2)-a*b*x(1)^2*x(2)-w0^2*r*x(1)^3+f*cos(omega*t)];
我遇到的问题就是每次w0 a b r omega设的值都不合适不出现分岔图。
望高手也试一试。如果试出来请告诉我呀。谢谢了。

[ 本帖最后由 咕噜噜 于 2009-8-18 08:19 编辑 ]

siyaoming

补充一下：w0最好设成1，omega最好设成接近1或者试稍微大于1别的随便了。

octopussheng

1 请检查程序中T的设置，应是有误的
2 取点也请考虑好，可与论坛中已有的程序对比一下

siyaoming

T怎么设置啊。改成T=2*pi/omega吗？
取点;tspan=0:T/25:1000;y(1000:2000:end,1)'是不是前面式中的T/25与
                                     后面式中的1000：2000：end
                                                                           25与2000不符合啊。

siyaoming

这是一时滞方程的分岔图。
clf
w0=1;a=-1;b=-0.3;r=0.2;omega=1.1;g1=1;g2=1;f=0.1:0.1:5;
lags=[pi/2,2*pi/3];T=2*pi/omega;
y1=[];G=[];y2=[];
for h=1:length(f)
sol=dde23(@yanchi,lags,[0.1 0],[0:T/100:40*T],f(h));
y1=[y1,sol.y(1,2000:100:end)];y2=[y2,sol.y(2,2000:100:end)];
m=length(y1);
for j=1:m/h
G=[G,f(h)];   
end;p=length(G);
N=length(f);
Over_to_step=N-i;
end
figure(1)
plot(G,y1,'k.','markersize',3)
hold on


function dy=yanchi(t,y,Z,f)
w0=1;
a=-1;
b=-0.3;
r=0.2;
omega=1.1;
g1=1;
g2=1;
f=0.1:0.1:5;
dy=[y(2);-w0^2*y(1)-a*y(2)-a*b*y(1)^2*y(2)-w0^2*r*y(1)^3+f*cos(omega*t)+g1*Z(1,1)+g2*Z(2,2)]

无水1324

另外，你也可以看看薛定宇的一本matlab书，里面也有一个很好的例子，我不知道适合你不？

siyaoming

楼上的院长啊，我去哪里看这本书呢？我没有。
帮帮我吧。

VibrationMaster

很恶心这种“跪求”。如果问道题目就这样，将来怎会有出息？

siyaoming

好的，下次不会了。盼望你给我一个解答。谢谢。

zhanghb

这是我的程序，可以参考一下
clear
x0=rand(1,4);
range=[0.1:0.2:10.1];
for tau=range;
sol=dde23('td',tau,x0,[0,2000]);
step=0.1;
x=0:step:2000;
y=deval(sol,x,1);
yp=deval(sol,x,2);
Y=[y;yp]';
tau
for k=1900/step:length(Y)
if abs(Y(k,2))<0.05 && Y(k,1)>0
    plot(tau,Y(k,1),'.','markersize',3)
    hold on;
end
end
end
xlabel('\it\tau')
ylabel('\itx_{1}')
title('\itBifurcation Diagram')