龙格库塔法画的分岔图，相图，poincare图对不上
第一个图是用四阶龙格库塔法画的分岔图的部分放大，omiga=1.288从图上看应该是单周期，第三个图是用扫频法做的poincare映射图，确实是只有一个点，第二个图是四阶龙格库塔做的相图，从数据对应看的确实是一个周期画出这么个东西，但是单周期对应的不应该就是个圈吗，这个看起来怎么像双倍周期对应的相图呢？求大佬解答。下面是程序（只是求解过程部分，输入未知量以及求解是用的另一个m文件）
分岔图程序：
function qiujie(omigab,omigae,v0,v1,nx,n,m,path1)
clc;
ad=cd;
omiga1=omigab:(omigae-omigab)/n:omigae;
omiga=omiga1*(2*pi);
options=odeset('RelTol',1e-7);
path=;
filename=['ox1',num2str(omigab),'_',num2str(omigae),'----v0.',num2str(v0),' & v1.',num2str(v1),' & nx.',num2str(nx),'.dat'];
filename2=['ox2',num2str(omigab),'_',num2str(omigae),'----v0.',num2str(v0),' & v1.',num2str(v1),' & nx.',num2str(nx),'.dat'];
mkdir(path);
fileID=fopen(,'w+');
fclose(fileID);
fileID=fopen(,'a+');
fileID2=fopen(,'w+');
fclose(fileID2);
fileID2=fopen(,'a+');
%%循环部分，每个循环解一个方程
for j=1:length(omiga)
t3=clock;
tt=2*pi/omiga(j);
[~,y]=ode45(@fangcheng,0:tt/100:800*tt,,options,v0,v1,nx,omiga(j));
fprintf(fileID,'%g\t',omiga1(j));
fprintf(fileID,'%g\t',y(35000-500+m:100:80001-500+m,1));
fprintf(fileID,'\n');
fprintf(fileID2,'%g\t',omiga1(j));
fprintf(fileID2,'%g\t',y(35000-500+m:100:80001-500+m,3));
fprintf(fileID2,'\n');
l=j/length(omiga);
t4=clock;
t=(length(omiga)-j)*etime(t4,t3)/60;
['进度：',num2str(l*100),'%| 预计剩时：',num2str(t)]
end
fclose(fileID);
fclose(fileID2);
相图与poincare映射图程序：
function qiujie(omiga,v0,v1,nx,path1)
clc;
ad=cd;
options=odeset('RelTol',1e-7);
path=;
mkdir(path);
fileID=fopen(,'w+');
fclose(fileID);
fileID=fopen(,'a+');
tt=2*pi/omiga;
[~,y]=ode45(@fangcheng,0:tt/100:800*tt,,options,v0,v1,nx,omiga);
a=;
=size(a);
for p=1:1:m
for q=1:1:n
if q==n
fprintf(fileID,'%g\n',a(p,q));
else
fprintf(fileID,'%g\t',a(p,q));
end
end
end
fclose(fileID);
clear fileID m n a p q;
fileID=fopen(,'w+');
fclose(fileID);
fileID=fopen(,'a+');
a=;
=size(a);
for p=1:1:m
for q=1:1:n
if q==n
fprintf(fileID,'%g\n',a(p,q));
else
fprintf(fileID,'%g\t',a(p,q));
end
end
end
fclose(fileID);
clear fileID m n a;
end
方程：
function y=fangcheng(t,x,v0,v1,nx,omiga)
y=[x(2);(-0.3-0.084*(v0+v1*sin(omiga*t)))*x(2)+...
(-0.22+0.021*(v0+v1*sin(omiga*t)))*x(4)+...
(-4.71525-0.0337*nx+0.2*(v0+v1*sin(omiga*t))+...
0.0523143*(v0+v1*sin(omiga*t))*(v0+v1*sin(omiga*t))-0.042*omiga*v1*cos(omiga*t))*x(1)+...
(8.13643-0.00475639*nx+1.0368434*(v0+v1*sin(omiga*t))+...
0.451728*(v0+v1*sin(omiga*t))*(v0+v1*sin(omiga*t))+0.01*omiga*v1*cos(omiga*t))*x(3)-...
936734*x(1)*x(3)*x(3)+...
245204*x(1)*x(1)*x(3)-...
76525.7*x(1)*x(1)*x(1)+...
69499.7*x(3)*x(3)*x(3);
x(4);(-0.2216-0.20235*(v0+v1*sin(omiga*t)))*x(2)+...
(-0.685-0.0978405*(v0+v1*sin(omiga*t)))*x(4)+...
(9.36946-0.0047564*nx-0.59457*(v0+v1*sin(omiga*t))+...
0.0814*(v0+v1*sin(omiga*t))*(v0+v1*sin(omiga*t))-0.1*omiga*v1*cos(omiga*t))*x(1)+...
(-181.568-0.255393*nx+0.2393*(v0+v1*sin(omiga*t))+...
0.836227*(v0+v1*sin(omiga*t))*(v0+v1*sin(omiga*t))-0.0489203*omiga*v1*cos(omiga*t))*x(3)+...
208499*x(1)*x(3)*x(3)-...
936734*x(1)*x(1)*x(3)+...
81734.6*x(1)*x(1)*x(1)-3863660*x(3)*x(3)*x(3)];
关于输入问题，分岔图输入的是hz，然后图的横坐标也是hz；但是相图和映射图输入的是弧度，就是1.288*2*pi，所以应该不是输入的问题
页:
[1]