请教各位高手一个受迫duffing方程的分叉图画法 急用谢谢

dekister

d(d(R)/dt)dt-(1/4)*v.^2*R-g*R.^3=V*R*(cos(t)).^2-r*v*d(R)/dt
V是控制参数（2.8~3）
v=2.03
g=-0.75
r=0.25
初始值[0.01，0]
请各位高手帮帮忙，小弟不胜感激谢谢！

gghhjj

给一系列的控制参数

然后针对每个控制参数用ode45求解方程的周期解

周期解收敛后，每个周期取一个点（和庞加莱图的去点方法一样）

然后将这一系列控制参数下所取的点组合在一起画图即可

程序还是自己写吧

gghhjj

下面给个例子

1. function dx=Duffing(t,x);
   
2. % Duffing方程
   
3. % dx=Duffing(t,[x;y])
   
4. %   t-时间，x,y-为自变量
   
5. %   eg: dx=Duffing(10,[0;0;39])
   
6. %
   
7. % 方程如下:
   
8. %   dx=y
   
9. %   dy=-r*y-x^3+b*cos(t)
   
10. %   r=0.3 b=39
    
11. %
    
12. % Author's email: ustb03-07@yahoo.com.cn
    
13. %
    
14. r=0.3;
    
15. b=x(3);
    
16. dx(1,1)=x(2);
    
17. dx(2,1)=-r*x(2)-x(1)^3+b*cos(t);
    
18. dx(3,1)=0;

复制代码

gghhjj

画相图

1. [T,Y]=ode45('Duffing',50,[0;0;39]);
   
2. plot(Y(:,1),Y(:,2));

复制代码

gghhjj

分岔图

1.    for p=linspace(6,13,280);
   
2.        [T,Y]=ode45('Duffing',10,[3;0;p]);
   
3.        [T,Y]=ode45('Duffing',100,Y(end,:));
   
4.        for k=2:length(Y)
   
5.           f=k-1;
   
6.           y=1;
   
7.           if Y(k,1)<0
   
8.              if Y(f,1)>0
   
9.                y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));
   
10.              end
    
11.           else
    
12.              if Y(f,1)<0
    
13.                y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));
    
14.              end
    
15.           end
    
16.           if y<0
    
17.              Z=[Z p+y*i];
    
18.           end
    
19.        end
    
20.    end

复制代码

mechanic05

前辈在吗？你给的分岔程序，中间部分看不懂，好象道理懂，但具体怎么实现还是没看懂．能否把中间部分给解释一下？我是初学者．我直接运行你的程序
for p=linspace(6,13,280);
       [T,Y]=ode45('Duffing',10,[3;0;p]);
       [T,Y]=ode45('Duffing',100,Y(end,:));
       for k=2:length(Y)
          f=k-1;
          y=1;
          if Y(k,1)<0
             if Y(f,1)>0
               y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));
             end
          else
             if Y(f,1)<0
               y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));
             end
          end
          if y<0
             Z=[Z p+y*i];
          end
       end
   end
plot(b,y)
也不通，请指教．非常感谢！！:@)

[ 本帖最后由 无水1324 于 2007-7-11 14:47 编辑 ]

pmsl

这个分岔图不太好吧？能看出分岔规律来吗？

无水1324

图不怎么好，但是程序有一定的借鉴意义

zx808

看看mark一下，为什么我还是不能看附件呢