软弹簧Duffing系统分岔，跪求高手指点迷津
软弹簧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=;
period=2*pi/w;
k=0;
YY1=[];
step=2*pi/100;%步长。
for fp=range
y0=;
fp
k=k+1;
%discard the first 60 periodic data;
%除去前面60个周期的数据，并将最后的结果作为下一次积分的初值
tspan=;
=ode45(@system_eq,tspan,y0);
y0=Y(end,:);
j=1;
for i=60:200
tspan=;
=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 发表于 2014-3-18 09:01
没细看，首先有个疑问，从主程序看你的方程应该是三个
但函数脚本中却只给出了两个，粘贴的时候漏了，还是 ...
还有一个，同样也是个软弹簧的 planet 发表于 2014-3-18 10:00
还有一个，同样也是个软弹簧的
建议补全代码，程序跑起来才容易找到问题 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=;
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);
=ode45(@system_eq,tspan,y0);
i=floor(n*400/M):N:n; %舍弃前40%的数据
hold on;
plot(fp,y(i,1),'r.','markersize',5);
%xlim();
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
总是不能出现我想要的结果，请各位老师帮忙看看，是程序有问题还是参数的问题
gghhjj 发表于 2014-3-18 09:01
没细看，首先有个疑问，从主程序看你的方程应该是三个
但函数脚本中却只给出了两个，粘贴的时候漏了，还是 ...
程序我贴出来了，老师帮忙看看 目前在学混沌，接触到Duffing系统模型，对弱信号识别监测很好。 猫头鹰先生 发表于 2014-4-3 19:47
目前在学混沌，接触到Duffing系统模型，对弱信号识别监测很好。
是的，软Duffing比较少 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个，你给了三个，略加调整后如下：
function dx=system_eq(t,x)
global wd fp k1 k2 k3 k4 k5 m xi w0
dx = zeros(2,1);
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
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.0;
xi=0.01;
w0=sqrt(k1/m);
a1=k2/m;
a2=k3/m;
a3=k4/m;
a4=k5/m;
wd=w0;
y0=;
N=500;
M=1000;
nn=21;
T=2*pi/wd;
%T=2;
tspan=0:T/N:M*T;
n=length(tspan);
for fp=linspace(0,1,nn)
display(fp);
=ode45(@system_eq,tspan,y0);
i=floor(n*400/M):N:n; %舍弃前40%的数据
hold on;
plot(fp,y(i,1),'r.','markersize',5);
%xlim();
end
运行结果并非你所说的直线，此处暂且不论你的结果是否正确，单纯从代码上看，求解思路应该是没有问题
至于结果是否正确需要结合你的模型特点进行分析判断
最好能够仔细看看个参数下的相轨迹，庞加莱图等 gghhjj 发表于 2014-4-8 15:00
至于结果是否正确需要结合你的模型特点进行分析判断
最好能够仔细看看个参数下的相轨迹，庞加莱图等
您的这个分岔图，跟我画得结果一样，不知道这是否是分岔图 planet 发表于 2014-4-8 22:48
您的这个分岔图，跟我画得结果一样，不知道这是否是分岔图
是分岔图，但至于分岔图本身是否正确还需要根据你的模型作进一步分析
看了一下0.5参数下的庞加莱图和相轨迹，你的该参数下，系统像进入的混沌状态 我想问一下 杜芬方程的幅频响应 你有研究吗 有的话 是怎么获得的 zzqmessi 发表于 2014-4-10 15:37
我想问一下 杜芬方程的幅频响应 你有研究吗 有的话 是怎么获得的
这个论坛有很多吧 planet 发表于 2014-4-11 09:21
这个论坛有很多吧
看了那么多 也没有看到对于跳跃现象解释的一个比较好的程序也没有找到 自己理论求解得到的结果也不吻合 就有点郁闷
页:
[1]
2