软弹簧Duffing系统分岔，跪求高手指点迷津

planet

软弹簧Duffing系统方程如下：function dx=system_eq(t,X)
dx=zeros(3,1);
dx(1)=y;

dx(2)=fp*sin(w*t)-x(1)+x(1)^3-2*xi*y;
分岔图程序如下：
global xi
range=[0.1:0.01:1];
period=2*pi/w;
k=0;
YY1=[];
step=2*pi/100;  %步长。
for fp=range
      y0=[0 0.001 0];
    fp
    k=k+1;
    %discard the first 60 periodic data;
    %除去前面60个周期的数据，并将最后的结果作为下一次积分的初值
    tspan=[0:step:60*period];
    [t,Y]=ode45(@system_eq,tspan,y0);
    y0=Y(end,:);
    j=1;
    for i=60:200
        tspan=[i*period:step:(i+1)*period];
        [t,Y]=ode45(@system_eq,tspan,y0);
        YY1(k,j)=Y(end,1);   % get the omega data from every period end
        j=j+1;               %取出每一个周期内的第一个解的最后一个值。
        y0=Y(end,:);
    end
end
bifdata=YY1(:,end-51:end);
plot(range,bifdata,'k.','markersize',1);


为什么我的分岔图是附件中的样子？参数已变化，忘老师们指教

gghhjj

没细看，首先有个疑问，从主程序看你的方程应该是三个
但函数脚本中却只给出了两个，粘贴的时候漏了，还是另有原因

planet

gghhjj 发表于 2014-3-18 09:01
没细看，首先有个疑问，从主程序看你的方程应该是三个
但函数脚本中却只给出了两个，粘贴的时候漏了，还是 ...

还有一个，同样也是个软弹簧的

gghhjj

planet 发表于 2014-3-18 10:00
还有一个，同样也是个软弹簧的

建议补全代码，程序跑起来才容易找到问题

planet

global wd fp k1 k2 k3 k4 k5 m xi a1 a2 a3 a4 w0
k1=2;
k2=10;
k3=-1;
k4=1;
k5=(w0^2+a3)*a2/(w0^2+a1);
m=1;
xi=0.01;
w0=sqrt(k1/m);
a1=k2/m;
a2=k3/m;
a3=k4/m;
a4=k5/m;
wd=w0;
y0=[0 0 0];
N=200;
M=1000;
nn=100;
T=2*pi/wd;
%T=2;
tspan=0:T/N:M*T;
n=length(tspan);
for fp=linspace(0,1,nn);
    display(fp);
    [t,y]=ode45(@system_eq,tspan,y0);
    i=floor(n*400/M):N:n;    %舍弃前40%的数据
    hold on;
    plot(fp,y(i,1),'r.','markersize',5);
    %xlim([0.9 1.2]);

end
这是我系统的方程
function dx=system_eq(t,X)
dx=zeros(2,1);
global wd fp k1 k2 k3 k4 k5 m xi w0
dx(1)=x(2);
if x(1)>0.0;
    dx(2)=fp*sin(wd*t)-k1/m*x(1)-k4/m*x(1)-k5/m*x(1)^3-2*xi*w0*x(2);
else
    dx(2)=fp*sin(wd*t)-k1/m*x(1)-k2/m*x(1)-k3/m*x(1)^3-2*xi*w0*x(2);
end
总是不能出现我想要的结果，请各位老师帮忙看看，是程序有问题还是参数的问题

planet

gghhjj 发表于 2014-3-18 09:01
没细看，首先有个疑问，从主程序看你的方程应该是三个
但函数脚本中却只给出了两个，粘贴的时候漏了，还是 ...

程序我贴出来了，老师帮忙看看

猫头鹰先生

目前在学混沌，接触到Duffing系统模型，对弱信号识别监测很好。

planet

猫头鹰先生 发表于 2014-4-3 19:47
目前在学混沌，接触到Duffing系统模型，对弱信号识别监测很好。

是的，软Duffing比较少

gghhjj

planet 发表于 2014-4-2 11:20
global wd fp k1 k2 k3 k4 k5 m xi a1 a2 a3 a4 w0
k1=2;
k2=10;

你的程序略有瑕疵，修改主要是函数变量不对，另外初值应该是2个，你给了三个，略加调整后如下：

1. function dx=system_eq(t,x)
   
2. global wd fp k1 k2 k3 k4 k5 m xi w0
   
3. dx = zeros(2,1);
   
4. dx(1) = x(2);
   
5. if x(1) > 0.0;
   
6.     dx(2) = fp * sin(wd * t) - k1/m * x(1) - k4 / m * x(1) - k5 / m * x(1) ^ 3 - 2 * xi * w0 * x(2);
   
7. else
   
8.     dx(2) = fp * sin(wd * t) - k1/m * x(1) - k2 / m * x(1) - k3 / m * x(1) ^ 3 - 2 * xi * w0 * x(2);
   
9. end

复制代码

1. global wd fp k1 k2 k3 k4 k5 m xi a1 a2 a3 a4 w0
   
2. k1=2;
   
3. k2=10;
   
4. k3=-1;
   
5. k4=1;
   
6. k5=(w0^2+a3)*a2/(w0^2+a1);
   
7. m=1.0;
   
8. xi=0.01;
   
9. w0=sqrt(k1/m);
   
10. a1=k2/m;
    
11. a2=k3/m;
    
12. a3=k4/m;
    
13. a4=k5/m;
    
14. wd=w0;
    
15. y0=[0 0];
    
16. N=500;
    
17. M=1000;
    
18. nn=21;
    
19. T=2*pi/wd;
    
20. %T=2;
    
21. tspan=0:T/N:M*T;
    
22. n=length(tspan);
    
23. 
    
24. for fp=linspace(0,1,nn)
    
25.     display(fp);
    
26.     [t,y]=ode45(@system_eq,tspan,y0);
    
27.     i=floor(n*400/M):N:n;    %舍弃前40%的数据
    
28.     hold on;
    
29.     plot(fp,y(i,1),'r.','markersize',5);
    
30.     %xlim([0.9 1.2]);
    
31. end

复制代码

运行结果并非你所说的直线，此处暂且不论你的结果是否正确，单纯从代码上看，求解思路应该是没有问题

gghhjj

至于结果是否正确需要结合你的模型特点进行分析判断
最好能够仔细看看个参数下的相轨迹，庞加莱图等

planet

gghhjj 发表于 2014-4-8 15:00
至于结果是否正确需要结合你的模型特点进行分析判断
最好能够仔细看看个参数下的相轨迹，庞加莱图等

您的这个分岔图，跟我画得结果一样，不知道这是否是分岔图

gghhjj

planet 发表于 2014-4-8 22:48
您的这个分岔图，跟我画得结果一样，不知道这是否是分岔图

是分岔图，但至于分岔图本身是否正确还需要根据你的模型作进一步分析
看了一下0.5参数下的庞加莱图和相轨迹，你的该参数下，系统像进入的混沌状态

zzqmessi

我想问一下 杜芬方程的幅频响应 你有研究吗 有的话 是怎么获得的

planet

zzqmessi 发表于 2014-4-10 15:37
我想问一下 杜芬方程的幅频响应 你有研究吗 有的话 是怎么获得的

这个论坛有很多吧

zzqmessi

planet 发表于 2014-4-11 09:21
这个论坛有很多吧

看了那么多 也没有看到对于跳跃现象解释的一个比较好的  程序也没有找到 自己理论求解得到的结果也不吻合 就有点郁闷