ruichard101 发表于 2008-4-6 11:39

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

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=
for h=1:length(gamma)
    T=2*pi/(w0*gamma(h));
    =ode45(@myfun_1,,x0,[],gamma(h));
    y(:,1)
    plot(gamma(h),y(18000:100:end,1),'k');hold on
end


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

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

无水1324 发表于 2008-4-6 11:49

回复 楼主 的帖子

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 发表于 2008-4-6 13:24

回复 2楼 的帖子

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

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

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

无水1324 发表于 2008-4-6 14:26

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)];



clc;
clear;
w0=25;
gamma=0.6;
x0=;
T=2*pi/(w0*gamma);
=ode45('myfun_1',,x0,[],gamma);
plot(y(1800:100:end,1),y(1800:100:end,2),'k');



现在应该可以运行了

ruichard101 发表于 2008-4-6 15:20

回复 4楼 的帖子

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

xinwilliam 发表于 2008-4-6 16:07

加个循环就行了

无水1324 发表于 2008-4-6 17:59

回复 4楼 的帖子

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

无水1324 发表于 2008-4-6 18:02

哎，........

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

octopussheng 发表于 2008-4-6 18:24

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

ruichard101 发表于 2008-4-6 19:15

回复 8楼 的帖子

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

octopussheng 发表于 2008-4-6 19:19

把错误提示贴一下吧！

ruichard101 发表于 2008-4-6 19:32

回复 11楼 的帖子

??? Input argument "gamma" is undefined.

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

这就是错误提示。

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

无水1324 发表于 2008-4-6 19:42

回复 12楼 的帖子

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

octopussheng 发表于 2008-4-6 19:49

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


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

ruichard101 发表于 2008-4-6 19:54

回复 13楼 的帖子

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

我改成‘myfun_1’以后能运行了。是不是得改x0的初值啊！初值是不是得取系统方程的一个奇点呢？
辛苦您了呵呵！:handshake
页: [1] 2
查看完整版本: 请高手们帮我看看这个碰磨方程的分叉程序