请高手们帮我看看这个碰磨方程的分叉程序

ruichard101

function dy=myfun_1(t,x)
global gamma;
q=3;
r=2;
beta=4;
epsilon=0.0025;
f=0.0025;
delta=0.16;
u=0.0845;
theta=0;
w0=25;
g=9.8;
dy=[x(2);
    -2*epsilon*x(2)-x(1)-beta*(sqrt(x(1)^2+x(3)^2)-delta)^(q/r)/sqrt(x(1)^2+x(3)^2)*(x(1)-f*x(3))+u*gamma^2*cos(gamma*w0*t+theta);
    x(4);
    -2*epsilon*x(4)-x(3)-beta*(sqrt(x(1)^2+x(3)^2)-delta)^(q/r)/sqrt(x(1)^2+x(3)^2)*(f*x(1)+x(3))+u*gamma^2*sin(gamma*w0*t+theta)-g/(w0^2)];

clc;
clear;
global gamma
w0=25;
gamma=0.6:0.01:1.3;
x0=[0;0;0;0]
for h=1:length(gamma)
    T=2*pi/(w0*gamma(h));
    [t y]=ode45(@myfun_1,[0:T/100:200*T],x0,[],gamma(h));
    y(:,1)
    plot(gamma(h),y(18000:100:end,1),'k');hold on
end


请各位大虾们帮我看看这个关于非线性碰磨的分叉图的程序，哪地方有问题！谢谢了！

[ 本帖最后由 ruichard101 于 2008-4-6 11:46 编辑 ]

无水1324

function dy=myfun_1(t,x,flag,gamma)
q=3;
r=2;
beta=4;
epsilon=0.0025;
f=0.0025;
delta=0.16;
u=0.0845;
theta=0;
w0=25;
g=9.8;
dy=[x(2);
    -2*epsilon*x(2)-x(1)-beta*(sqrt(x(1)^2+x(3)^2)-delta)^(q/r)/sqrt(x(1)^2+x(3)^2)*(x(1)-f*x(3))+u*gamma^2*cos(gamma*w0*t+theta);
    x(4);
    -2*epsilon*x(4)-x(3)-beta*(sqrt(x(1)^2+x(3)^2)-delta)^(q/r)/sqrt(x(1)^2+x(3)^2)*(f*x(1)+x(3))+u*gamma^2*sin(gamma*w0*t+theta)-g/(w0^2)];

ruichard101

无水您好！我按照你说的那样改了，运行出现了下面的错误：
??? Input argument "gamma" is undefined.

Error in ==> myfun_1 at 13
dy=[x(2);

我想问一下：初值是不是得取方程的一个奇点呢？这个方程我解了好久，就是得不到结果。我的主程序是不是有问题？我的原意是：积分步长是T/100，舍去前180个周期，只取后20个周期的解

无水1324

1. function dy=myfun_1(t,x,flag,gamma)
   
2. q=3;
   
3. r=2;
   
4. beta=4;
   
5. epsilon=0.0025;
   
6. f=0.0025;
   
7. delta=0.16;
   
8. u=0.0845;
   
9. theta=0;
   
10. w0=25;
    
11. g=9.8;
    
12. dy=[x(2);
    
13. -2*epsilon*x(2)-x(1)-beta*(sqrt(x(1)^2+x(3)^2)-delta)^(q/r)/sqrt(x(1)^2+x(3)^2)*(x(1)-f*x(3))+u*gamma^2*cos(gamma*w0*t+theta);
    
14. x(4);
    
15. -2*epsilon*x(4)-x(3)-beta*(sqrt(x(1)^2+x(3)^2)-delta)^(q/r)/sqrt(x(1)^2+x(3)^2)*(f*x(1)+x(3))+u*gamma^2*sin(gamma*w0*t+theta)-g/(w0^2)];

复制代码

1. 
   
2. clc;
   
3. clear;
   
4. w0=25;
   
5. gamma=0.6;
   
6. x0=[0.5;0.5;0.5;0.5];
   
7. T=2*pi/(w0*gamma);
   
8. [t y]=ode45('myfun_1',[0:T/100:200*T],x0,[],gamma);
   
9. plot(y(1800:100:end,1),y(1800:100:end,2),'k');
   

复制代码

现在应该可以运行了

ruichard101

无水大哥不好意思呵呵！我没有说清楚。我想做的是频率比（gamma）在0。6----1。3之间的响应的分叉图。gamma不是一固定值。请教一下应该怎么编写主程序呢？

xinwilliam

加个循环就行了

无水1324

我的意思是上面的程序可以运行了，因为我不可能跟你试那个分叉的程序。
你加一个循环就可以了

无水1324

哎，........

clc;
clear;
w0=25;
gamma=0.6:0.01:1;
x0=[0.5;0.5;0.5;0.5];
for h=1:length(gamma)
T=2*pi/(w0*gamma(h));
[t y]=ode45('myfun_1',[0:T/100:200*T],x0,[],gamma(h));
plot(gamma(h),y(1800:100:end,2),'k.');hold on
end

octopussheng

呵呵，无水，你辛苦了哦！

ruichard101

无水大哥，多次麻烦您实在不好意思:@$ 。我知道得加一个循环。我上面的主程序和您给的基本一样啦！~~可就是运行的时候老是出现我在三楼给出的错误啊！折腾好几天了，改来改去总是不行:@L
也谢谢xin和oct的指点！

octopussheng

把错误提示贴一下吧！

ruichard101

??? Input argument "gamma" is undefined.

Error in ==> myfun_1 at 12
dy=[x(2);

这就是错误提示。

gamma是个控制参数，怎么就没定义呀？

无水1324

我运行了没有出现你说的错误，但就是解出现了复数

octopussheng

无水，我看他还是没有把matlab求解微分方程的方法看懂！呵呵！


回复ruichard101：无水在4楼给你贴了两段程序，第一个程序是定义微分方程的程序，是不能直接运行的，你上面之所以会有问题，就是直接运行这个程序造成的；第二个程序才是计算求解的代码，应该运行第二个程序才能得到结果！

ruichard101

无水大哥，我仔细比较了一下。。。
@myfun_1和‘myfun_1’这两种表示不是都行吗？我用的是@myfun_1，所以导致这几天编的一个又一个程序就这么被自己“扼杀”了。。。:@(

我改成‘myfun_1’以后能运行了。是不是得改x0的初值啊！初值是不是得取系统方程的一个奇点呢？
辛苦您了呵呵！:handshake