庞加莱截面的c语言编程,画Lorenz庞加莱截面图的程序

% 画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('d\y','fontsize',14);

%

作此截面图时，其本质是当非线性系统进行角动量作用变换后，可以在环面上讨论系统的性质。

%

环面是指类似轮胎内胎的一个东西，而轨线在环面上运动，这种运动是复杂的，

%

包含绕环中心的运动(公转)和绕截面圆心的运动(类似自转)；如果这些有两种频率，

% 分别为a

,b，当a/b为有理数时，它们最终会首尾结合到一起，从而形成一个圈。

%

因此，截面上只会留下一个点，这个点是轨线穿越截面时留下的；当两种频率之比(旋转数)为无理数时，

%

两轨线永不相交，这时截面上留下一个圆的痕迹，所以所拟周期的Poincare截面图是一个圆。