马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
分岔图(部分放大)
相图(omiga=1.288)
poincare映射图(omiga=1.288)
第一个图是用四阶龙格库塔法画的分岔图的部分放大,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=[path1,'\ox\v0.',num2str(v0),' & v1.',num2str(v1),' & nx.',num2str(nx)];
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([ad,'\',path,'\',filename],'w+');
fclose(fileID);
fileID=fopen([ad,'\',path,'\',filename],'a+');
fileID2=fopen([ad,'\',path,'\',filename2],'w+');
fclose(fileID2);
fileID2=fopen([ad,'\',path,'\',filename2],'a+');
%%循环部分,每个循环解一个方程
for j=1:length(omiga)
t3=clock;
tt=2*pi/omiga(j);
[~,y]=ode45(@fangcheng,0:tt/100:800*tt,[0.001,0.001,0.001,0.001],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=[path1,' 1bi3 neigongzhen','\omiga.',num2str(omiga),' & v0.',num2str(v0),' & v1.',num2str(v1)];
mkdir(path);
fileID=fopen([ad,'\',path,'\','V1相图','.dat'],'w+');
fclose(fileID);
fileID=fopen([ad,'\',path,'\','V1相图','.dat'],'a+');
tt=2*pi/omiga;
[~,y]=ode45(@fangcheng,0:tt/100:800*tt,[0.001,0.001,0.001,0.001],options,v0,v1,nx,omiga);
a=[y(35000:80001,1),y(35000:80001,2)];
[m,n]=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([ad,'\',path,'\','V1映射图','.dat'],'w+');
fclose(fileID);
fileID=fopen([ad,'\',path,'\','V1映射图','.dat'],'a+');
a=[y(35000:100:80001,1),y(35000:100:80001,2)];
[m,n]=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)];
|