学习混沌: 庞加莱截面和倍周期分岔

以前的关于微分方程组的庞加莱截面和倍周期分岔程序有误，现在更正如下：


Z=[];
[T,Y]=ode45('Lorenz',[0,1000],[1;1;1;16;4;49.52]);
Y(:,3)=Y(:,3)-30;
for k=1:(length(Y)-1)
  if Y(k,3)<0
    if Y(k+1,3)>0
      % 用直线连接位于截面两侧的相邻离散点，近似得到穿过截面的点
       x=Y(k,1)-Y(k,3)*(Y(k+1,1)-Y(k,1))/(Y(k+1,3)-Y(k,3));
       y=Y(k,2)-Y(k,3)*(Y(k+1,2)-Y(k,2))/(Y(k+1,3)-Y(k,3));
       Z=[Z x+i*y];
    end

  % 如果只画从一侧穿过截面的图象，以下的if-end语句可删去
  else
    if Y(k+1,3)<0
       x=Y(k,1)-Y(k,3)*(Y(k+1,1)-Y(k,1))/(Y(k+1,3)-Y(k,3));
       y=Y(k,2)-Y(k,3)*(Y(k+1,2)-Y(k,2))/(Y(k+1,3)-Y(k,3));
       Z=[Z x+i*y];
     end

  end
end
plot(Z,'.','markersize',2)
title('Lorenz 系统的 Poincare 映像 z=30')
xlabel('x'),ylabel('y')


Z=[];
[T,Y]=ode45('Chen',[0,1000],[1;1;1;35;3;28]);
for k=1:(length(Y)-1)
  if Y(k,1)<0
    if Y(k+1,1)>0
       y=Y(k,2)-Y(k,1)*(Y(k+1,2)-Y(k,2))/(Y(k+1,1)-Y(k,1));
       z=Y(k,3)-Y(k,1)*(Y(k+1,3)-Y(k,3))/(Y(k+1,1)-Y(k,1));
       Z=[Z y+i*z];
    end
  else
    if Y(k+1,1)<0
       y=Y(k,2)-Y(k,1)*(Y(k+1,2)-Y(k,2))/(Y(k+1,1)-Y(k,1));
       z=Y(k,3)-Y(k,1)*(Y(k+1,3)-Y(k,3))/(Y(k+1,1)-Y(k,1));
       Z=[Z y+i*z];
     end
  end
end
plot(Z,'.','markersize',2)
title('Chen 系统的 Poincare 映像 x=0')
xlabel('y'),ylabel('z')


Z=[];
[T,Y]=ode45('Lorenz3',[0,5000],[0.5;0.5;2;.9;.1;2;0]);
Y(:,2)=Y(:,2)-.2;
for k=1:(length(Y)-1)
  if Y(k,2)<0
    if Y(k+1,2)>0
       x=Y(k,1)-Y(k,2)*(Y(k+1,1)-Y(k,1))/(Y(k+1,2)-Y(k,2));
       z=Y(k,3)-Y(k,2)*(Y(k+1,3)-Y(k,3))/(Y(k+1,2)-Y(k,2));
       Z=[Z x+i*z];
    end
  else
    if Y(k+1,2)<0
       x=Y(k,1)-Y(k,2)*(Y(k+1,1)-Y(k,1))/(Y(k+1,2)-Y(k,2));
       z=Y(k,3)-Y(k,2)*(Y(k+1,3)-Y(k,3))/(Y(k+1,2)-Y(k,2));
       Z=[Z x+i*z];
     end
  end
end
plot(Z,'.','markersize',2)
title('模拟 Lorenz 系统的 Poincare 映像 y=0.2')
xlabel('x'),ylabel('z')


Z=[]; 
for r=linspace(1,500,250);
  % 舍弃前面迭带的结果，用后面的结果画图
  [T,Y]=ode45('Lorenz',1,[1;1;1;10;8/3;r]);
  [T,Y]=ode45('Lorenz',50,Y(length(Y),:));
  Y(:,1)=Y(:,2)-Y(:,1);
  for k=2:length(Y)
    f=k-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));
         Z=[Z r+abs(y)*i];
       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));
         Z=[Z r+abs(y)*i];
       end
     end
   end
end 
plot(Z,'.','markersize',1) 
title('Lorenz映射分岔图')
xlabel('r'),ylabel('|y| where x=y')

学习混沌

博客归档

Saturday, October 21, 2006

庞加莱截面和倍周期分岔

No comments: