% 画Lorenz庞加莱截面图的程序
LorenzDiEqn2=inline(['[-10*(x(1)-x(2));',...
'-x(1)*x(3)+28*x(1)-x(2);',...
'x(1)*x(2)-8/3*x(3)]'],'t','x');
clc;
% function dx = LorenzDifEqn2(t,x);
% Lorenz系统微分方程
% 方程如下:
% dx = -σ*(x-y)
% dy = -x*z+r*x-y
% dz = x*y-b*z
% sigma = 10;
% b = 8/3;
% r = 20;
% d_x = -sigma*(x(1)-x(2));
% dy = -x(1)*x(3)+r*x(1)-x(2);
% dz = x(1)*x(2)-b*x(3);
% dx = [d_x;dy;dz];
% Poincare_section[绘制庞加莱截面图]
[t,x]=ode45(LorenzDiEqn2,[0,100],[0.1,0.1,10])
figure
plot3(x(:,1),x(:,2),x(:,3))
%figure
%plot(t,x(:,2))
%figure
%plot(t,x(:,3))
z0=28; % 选择z0=28这个截面
j = 0;
for k = 1:length(x(:,3))-1
d1 = x(k,3)-z0;
d2 = x(k+1,3)-z0;
if abs(d1)<1e-8
j = j+1;
X1(j) = x(k,1);
X2(j) = x(k,2);
continue;
end
if sign(d1)*sign(d2)<0
j = j+1;
Q=polyfit([x(k,3),x(k+1,3)],[x(k,1),x(k+1,1)],1);
X1(j)=polyval(Q,z0);
Q=polyfit([x(k,3),x(k+1,3)],[x(k,2),x(k+1,2)],1);
X2(j)=polyval(Q,z0);
end
end
figure
plot(X1,X2,'.');
xlabel('x','fontsize',14);
ylabel('dy','fontsize',14);
%
作此截面图时,其本质是当非线性系统进行角动量作用变换后,可以在环面上讨论系统的性质。
%
环面是指类似轮胎内胎的一个东西,而轨线在环面上运动,这种运动是复杂的,
%
包含绕环中心的运动(公转)和绕截面圆心的运动(类似自转);如果这些有两种频率,
% 分别为a
,b,当a/b为有理数时,它们最终会首尾结合到一起,从而形成一个圈。
%
因此,截面上只会留下一个点,这个点是轨线穿越截面时留下的;当两种频率之比(旋转数)为无理数时,
%
两轨线永不相交,这时截面上留下一个圆的痕迹,所以所拟周期的Poincare截面图是一个圆。