MATLAB混沌系统仿真：Lorenz系统和Rossler系统
MATLAB混沌系统仿真：Lorenz系统和Rossler系统
本文将介绍如何使用MATLAB进行混沌系统的仿真，以Lorenz系统和Rossler系统为例。通过龙格库塔四阶算法求解微分方程，并计算最大李亚普诺夫指数来判断系统的混沌特性。
一、混沌产生系统：Lorenz和Rossler
Lorenz系统是一个按Rayleigh-Bénard配置的大气对流简化模型，包含3个微分方程。Lorenz在1963年证明了这个系统中存在混沌。变量X、Y、Z分别表示循环流体的流速、上升和下降流体的温差和垂直温度剖面的畸变。其微分形式如下：
$$
\begin{align*}
\frac{dx}{dt} &= -\sigma(x-y) \
\frac{dy}{dt} &= rx - y - xz \
\frac{dz}{dt} &= -\beta z + xy
\end{align*}
$$
其中控制参数 $\sigma = 10$、$r = 28$、$\beta = \frac{8}{3}$ 时也叫经典Lorenz系统，本次展示也是用这些参数。经典参数下Lorenz系统最大李亚普诺夫指数计算得2.12，系统最终将走向混沌。同时洛伦兹方程是耗散的，所有的轨迹最终进入吸收域：
$$
\Omega = \left{ x \in \mathbb{R}^3 : rx^2 + \sigma y^2 + \sigma(z-2r)^2 < \frac{\beta^2 r^2}{\beta-1} \right}
$$
这就是那个蝴蝶型的区域，Lorenz吸引子的每一条轨迹在空间上都不重合，并且总不会超出这个区域。
Rossler系统是与人为构造的与Lorenz类似可以产生混沌的系统，并没有实际的物理意义，因此不做过多赘述。其方程如下：
$$
\begin{align*}
\frac{dx}{dt} &= -(y+z) \
\frac{dy}{dt} &= x + ay \
\frac{dz}{dt} &= b + z(x-c)
\end{align*}
$$
二、龙格库塔四阶算法
求解上述两个方程组，MATLAB提供ode函数求解，在精度要求不高时，与龙格库塔四阶算法差别不大。龙格库塔四阶算法在每个步长 $h$ 上使用 $k_1$、$k_2$、$k_3$、$k_4$ 四个分别代表开始斜率、考虑 $k_1$ 的中点斜率、考虑 $k_2$ 的中点斜率、考虑 $k_3$ 的末端斜率，取加权平均：
$$
y(n+1) = y(n) + \frac{h}{6} \cdot (k_1 + 2 \cdot k_2 + 2 \cdot k_3 + k_4)
$$
这个方法的求解具有4阶精度，跟DDE方法(实质为三阶变步长的龙格库塔算法)精度相近，但是运算更快。其代码如下：
function [t,y]=rk4(f,tspan,y0,n)
%四阶Runge-Kutta算法
%f:dy/dt=f(t,y);求解区间tspan:[t0,tf];初值y0=y(t0)=[x0,y0,z0];n采样点数
h=(tspan(2)-tspan(1))/n;
t=tspan(1):h:tspan(2);
y=zeros(length(y0),length(t));
y(:,1)=y0;
for m=2:length(t)
k1=f(t(m-1),y(:,m-1));
k2=f(t(m-1)+h/2,y(:,m-1)+k1*h/2);
k3=f(t(m-1)+h/2,y(:,m-1)+k2*h/2);
k4=f(t(m-1)+h,y(:,m-1)+k3*h);
y(:,m)=y(:,m-1)+(k1+2*k2+2*k3+k4)*h/6;
end
t=t';
y=y';
end
其中需要首先传入一个函数的句柄，函数可由如下形式定义：
function dy_dt=lorenz(t,y)
%y=[x,y,z]储存三个变量的矩阵
p=10;
r=28;
b=8/3;
dx_dt=p*(y(2)-y(1));
dy_dt=r*y(1)-y(2)-y(1)*y(3);
dz_dt=y(1)*y(2)-b*y(3);
dy_dt=[dx_dt;dy_dt;dz_dt];
end
Rossler函数的定义相同，从Lorenz形式的方程和代码很容易写出，这里不做累述。值得一提的是，如果不是做半导体激光器速率方程那样 时滞 的仿真，实际上只需要简单一行代码调用ode45函数就可以出结果，但是ode函数是 变步长 的，因此传入的需要传入固定距离的长矩阵(见下文)。
三、混沌吸引子轨迹
运行后(完整代码见附录)即可得到蝴蝶型的吸引子。使用常规的frame方法保存图像如下：(用view函数将图像稍稍转了下)。与自带的ODE函数运行的结果几乎一致。
Rossler吸引子如下：
四、最大李亚普诺夫指数
如何将chaos与“噪声”和“周期特别长的序列”中区别开，最大李指数是一种不错的方法。如果一个系统具有一个正的李指数，那么它就是混沌的；正的数值越大，系统越混沌。因此也有人(吕金虎，2002，如果没记错的话)认为 李指数的倒数 是 系统可预测时间的测度 ，超过这个值，系统将不可预测。
此处参考Wolf博士的程序(at the University of Texas at Austin)。他用BASGEN生成数据库，与用FET生成的时间序列一起应用于监测邻近轨迹的散度来估计最大李指数。FET和BASGEN实际上都使用了延时重建和相空间重构的方法，因此实际上对结果影响较大的参数也就主要是嵌入维度、延时、演化时间，此处分别设为3,10,20。对采样出16384个点的Lorenz系统进行估计。对轨道的检测结果如下如下：
可以看出，在大部分的演化时间里，两个轨道紧紧贴合，发散程度低。下面是FET的输出，第5列是李雅普诺夫指数的估计。第5列是轨道角修正数值。
程序预测出的李指数结果为2.15，与计算值接近，因此Lorenz系统是混沌的。
（2023.7更新）
当然，MATLAB自带的程序也可以实现重构相空间与计算最大李指数。使用MATLAB自带的Lorenz数据，频率fs为100。使用phaseSpaceReconstruction函数使用平均互信息重构数据的相空间，得到维度为3，重构延时为10。数据量较大时重构的速度慢。此处lyapunovExponent的算法同样可以表示闭合相空间的轨迹分离率，计算得1.62。李指数是有单位的，由频率 $f_s$ 决定，此处为 $s^{-1}$ 。如果两个点处于同一起点，那么经过平均1.62s后它们的轨道将完全分离。
%使用matlab自带的Lorenz数据
load('lorenzAttractorExampleData.mat','data','fs');
xdata = data(:,1);
[~,lag,dim] = phaseSpaceReconstruction(xdata)
eg=[80 150]
lyp=lyapunovExponent(xdata,fs,lag,dim,'ExpansionRange',eg);%取值
lyapunovExponent(xdata,fs,lag,dim,'ExpansionRange',eg);%绘图
参数ExpansionRange的选取，线性拟合越好的区间越合适。下图是一个特别大的区间[0 500]，无法线性拟合，计算结果偏差就很大。
附：完整程序
y0=[1;0;0];
tspan=[0,1e2];
[t1,y1]=rk4(@lorenz,tspan,y0,1e5);
[t2,y2]=ode45(@lorenz,tspan,y0);
plot3(y1(:,1),y1(:,2),y1(:,3));
figure();plot3(y2(:,1),y2(:,2),y2(:,3));
function [t,y]=rk4(f,tspan,y0,n)
%四阶Runge-Kutta算法
%f:dy/dt=f(t,y);求解区间tspan:[t0,tf];初值y0=y(t0)=[x0,y0,z0];n采样点数
h=(tspan(2)-tspan(1))/n;
t=tspan(1):h:tspan(2);
y=zeros(length(y0),length(t));
y(:,1)=y0;
for m=2:length(t)
k1=f(t(m-1),y(:,m-1));
k2=f(t(m-1)+h/2,y(:,m-1)+k1*h/2);
k3=f(t(m-1)+h/2,y(:,m-1)+k2*h/2);
k4=f(t(m-1)+h,y(:,m-1)+k3*h);
y(:,m)=y(:,m-1)+(k1+2*k2+2*k3+k4)*h/6;
end
t=t';
y=y';
end
function dy_dt=lorenz(t,y)
%y=[x,y,z]储存三个变量的矩阵
p=10;
r=28;
b=8/3;
dx_dt=p*(y(2)-y(1));
dy_dt=r*y(1)-y(2)-y(1)*y(3);
dz_dt=y(1)*y(2)-b*y(3);
dy_dt=[dx_dt;dy_dt;dz_dt];
end