非线性可视化（3）混沌系统

承接上一篇二维相图。

如果二维相平面中出现了交叉的轨线，则说明这个系统的维度很可能大于二维。

下面就以几个经典的系统作为示范。本章不涉及太多知识点，以展示为主。主要介绍三个经典的非线性混沌系统。

• 1 Lorenz系统

Lorenz系统是气象学家洛伦兹发现并提出的一个非线性系统，也是混沌学科的开端。在模拟大气流动时，洛伦兹发现初始的一个小小的误差，都会导致系统未来极大的变化。这种思想在20世纪60年代，给了那些物理学界中决定论者沉重的打击。洛伦兹也将这种不确定性，总结为“蝴蝶效应”。

这个系统可以被写为：

一般系统a=10，b=8/3，变化r值来观察系统的不同样子。下图分别展示了取不同r值所对应的xy平面的二维相轨线图：

其中r=20时，对应系统收敛到定点。r=28对应混沌。r=99.36对应倍周期。

因为二维系统在相平面上不会出现交叉，所以混沌、倍周期等现象都是在三维或者更高维才出现。正如前文所说的，对于混沌来说，三是个神奇的数字。

它们在三维空间中的轨迹图为：

中间的那个图就是经典的洛伦兹吸引子图。

• 2 Rossler系统

Rossler系统是Rössler本人在70年代提出的一个非线性系统，和前面的Lorenz系统相比更为简单，但是却依然拥有复杂的非线性行为。

它可以写作：

下图绘制了a=0.1，b=0.1，改变不同的c绘制的轨迹图。

其中周期2指的是每两个波形一个循环，系统转2圈回到同一个点。随着c的增大，系统由周期1到了周期2，之后突然增加大周期4，再之后以越来越快的速度增大到周期8甚至更高，最后密密麻麻的周期似乎再也不会循环，变成了混沌。

这种周期越来越多逐渐变为混沌的现象，叫做倍周期分岔现象。是系统由有序变为无序混沌常见的一种方式。

• 3 duffing方程

duffing方程也是以 Georg Duffing命名的一个非线性方程。它是基于强迫振动的单摆所提出的方程，它提出的时间非常早，但是被拿来做混沌研究还是比较晚的。由于它背后有着非常明显与简单的物理模型，所以甚至可以做实验去观察这个方程的非线性[3]。方程的形式为：

与前面两个方程不同，duffing方程存在一个强迫振动项，带有时间t，所以不属于自治系统。可以看到虽然系统是二阶的，但仍然具有非常复杂的非线性。

如下图，固定激励的振幅频率r和w，改变阻尼d。

可以看到随着阻尼d的增大，系统由混沌变为2周期，又变为了单周期运动。

对于这种一团乱麻的混沌现象，只观察轨迹图并不能看到什么规律。下一章节我们将引入一种新的观测方法——庞佳莱截面法。

后面附上代码：

clc
clear
close all
%% 洛伦兹吸引子
h=1e-3;
x0=0:h:40;
[y1,~]=ODE_RK4_hyh(x0,h,[1;4;20],{'Lorenz',[10,8/3,20]});
Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);


figure(1)
subplot(1,3,1)
plot(Lx,Ly);title('r=20')
figure(2)
subplot(1,3,1)
plot3(Lx,Ly,Lz);view([51,30]);title('r=20')


[y1,~]=ODE_RK4_hyh(x0,h,[-13;-2;41],{'Lorenz',[10,8/3,28]});
Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);


figure(1)
subplot(1,3,2)
plot(Lx,Ly);title('r=28')
figure(2)
subplot(1,3,2)
plot3(Lx,Ly,Lz);view([51,30]);title('r=28')


[y1,~]=ODE_RK4_hyh(x0,h,[1;4;67],{'Lorenz',[10,8/3,99.36]});
Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);


figure(1)
subplot(1,3,3)
plot(Lx,Ly);title('r=99.36')
figure(2)
subplot(1,3,3)
plot3(Lx,Ly,Lz);view([51,30]);title('r=99.36')


%% Rossler吸引子
h=2e-3;
x0=0:h:180;
[y1,~]=ODE_RK4_hyh(x0,h,[4.9;-5;0.07],{'Rossler',[0.1,0.1,4]});
Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);


figure(4)
subplot(2,2,1)
plot3(Lx,Ly,Lz);view([51,30]);title('c=4 周期1')


[y1,~]=ODE_RK4_hyh(x0,h,[9.1;-5;0.17],{'Rossler',[0.1,0.1,6]});
Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);


figure(4)
subplot(2,2,2)
plot3(Lx,Ly,Lz);view([51,30]);title('c=6 周期2')


[y1,~]=ODE_RK4_hyh(x0,h,[12.8;-5;0.277],{'Rossler',[0.1,0.1,8.5]});
Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);


figure(4)
subplot(2,2,3)
plot3(Lx,Ly,Lz);view([51,30]);title('c=8.5 周期4')


[y1,~]=ODE_RK4_hyh(x0,h,[12.8;-5;0.277],{'Rossler',[0.1,0.1,9]});
Lx=y1(1,:);Ly=y1(2,:);Lz=y1(3,:);


figure(4)
subplot(2,2,4)
plot3(Lx,Ly,Lz);view([51,30]);title('c=9 混沌')


%% Duffing吸引子
h=2e-3;
x0=0:h:180;
[y1,~]=ODE_RK4_hyh(x0,h,[1;0.5],{'Duffing',[1.15,1,1]});
Lx=y1(1,:);Ly=y1(2,:);
figure(6)
subplot(1,3,1)
plot(Lx,Ly);title('d=1.15')


[y1,~]=ODE_RK4_hyh(x0,h,[0.8;0.75],{'Duffing',[1.35,1,1]});
Lx=y1(1,:);Ly=y1(2,:);
figure(6)
subplot(1,3,2)
plot(Lx,Ly);title('d=1.35')


[y1,~]=ODE_RK4_hyh(x0,h,[0.7;0.73],{'Duffing',[1.5,1,1]});%[0.7;0.73]
Lx=y1(1,:);Ly=y1(2,:);
figure(6)
subplot(1,3,3)
plot(Lx,Ly);title('d=1.5')


function [F,Output]=Fdydx(x,y,Input)
%形式为Y'=F(x,Y)的方程，参见数值分析求解常系数微分方程相关知识
%高次用列向量表示，F=[dy(1);dy(2)];y(1)为函数，y(2)为函数导数
switch Input{1}
    case 'Lorenz'
        a=Input{2}(1);b=Input{2}(2);r=Input{2}(3);
        dy(1)=a*(y(2)-y(1));
        dy(2)=r*y(1)-y(2)-y(1)*y(3);
        dy(3)=y(1)*y(2)-b*y(3);
        F=[dy(1);dy(2);dy(3)];
    case 'Rossler'
        a=Input{2}(1);b=Input{2}(2);c=Input{2}(3);
        dy(1)=-y(2)-y(3);
        dy(2)=y(1)+a*y(2);
        dy(3)=b+y(3)*(y(1)-c);
        F=[dy(1);dy(2);dy(3)];
    case 'Duffing'
        d=Input{2}(1);r=Input{2}(2);w=Input{2}(3);
        dy(1)=y(2);
        dy(2)=-y(1)^3+y(1)-d*y(2)+r*cos(w*x);
        F=[dy(1);dy(2)];
end
Output=[];
end


function [y,Output]=ODE_RK4_hyh(x,h,y0,Input)
%4阶RK方法
%h间隔为常数的算法
y=zeros(size(y0,1),size(x,2));
y(:,1)=y0;
for ii=1:length(x)-1
    yn=y(:,ii);
    xn=x(ii);
    [K1,~]=Fdydx(xn    ,yn       ,Input);
    [K2,~]=Fdydx(xn+h/2,yn+h/2*K1,Input);
    [K3,~]=Fdydx(xn+h/2,yn+h/2*K2,Input);
    [K4,~]=Fdydx(xn+h  ,yn+h*K3  ,Input);
    y(:,ii+1)=yn+h/6*(K1+2*K2+2*K3+K4);
end
Output=[];
end

参考资料：

[1] 微分方程、动力系统与混沌导论[M]

[2] Duffing equation Wiki 

[3] 计算物理基础-第10章第77讲（北京师范大学）（中国大学MOOC）计算物理基础_北京师范大学_中国大学MOOC(慕课) (icourse163.org)