基于混沌系统的文本加密算法研究（二）——经典混沌映射

---

前言

上一篇文章介绍了混沌以及混沌加密的一些基础知识，本文将介绍经典的混沌映射。
（传送门：基于混沌系统的文本加密算法研究（一）——混沌及混沌加密的基础知识）

---

本章将介绍经典的混沌映射，包括Logistic映射、Henon映射以及Lorenz映射，通过数值仿真，对混沌系统的特征及重要参数进行研究。（划重点，有完整的代码哦~~）在此基础上，后面将使用经典混沌映射对文本进行加密。

（所有的代码均在文章的最后给出）

一、一维Logistic混沌映射

Logistic映射是一个经典的一维混沌映射。该映射的方程如下：$$x_{n+1}=\mu x_n(1-x_n)\text{\hspace{1em}(1)}$$ 其中，$\mu$为控制参数，满足$\mu \in \left(0,4\right]，x\in \left(0.1\right)$。
当$\mu \in \left(0,4\right]$时，一维Logistic混沌映射的分岔图如下图所示。

由图3-1可知，当$0\le \mu \le 1$时，系统仅有一个不动点$x_0=0$，呈现静止状态；当$1<\mu \le 3$时，对于每个确定的$\mu$，系统只有一个非0的稳定状态，呈现周期1解；当$3<\mu \le 3.40$时，系统呈现周期2解；当$3.4<\mu \le 3.54$时，系统呈现周期4解。当控制参数$\mu$持续增加时，Logistic系统将依次呈现周期8解、周期16解以及周期32解等。即随着参数$\mu$的增加，系统会发生倍周期分岔的现象。当参数$\mu$超过极限值${\mu }_{\infty }=3.57\cdots \cdots$时，系统将进入混沌区，出现混沌现象，x的取值变得随机。但混沌区域中也存在一些空白带，称为“窗口”，对于窗口中的$\mu$，系统处于非混沌状态。即并非所有大于${\mu }_{\infty }$的映射都出现混沌，如当$\mu =3.836$时，系统将呈现周期3解。
进一步考察当$\mu$变化时Logistic映射的时间序列变化情况，分别取$\mu$为1、3.25、3.5和3.6，得到x的时间序列如下图所示。

由图3-2可知，随着$\mu$的增加，Logistic系统通过倍周期分岔进入混沌状态。
对$\mu \in \left(0,4\right]$，计算得到Logistic映射的Lyapunov特征指数谱如下图所示。

从图3-3可以看出，一维Logistic混沌映射的Lyapunov指数在参数$\mu >3.57$后大于0，即系统具有正的Lyapunov特征指数，系统进入混沌状态。

二、二维Henon混沌映射

Henon映射是一个经典的二维离散混沌映射。其方程如下：$$\{\begin{array}{l}{x}_{n+1}=1+y_n-a{x}_n^2\\ y_{n+1}=b{x}_n\end{array}\text{\hspace{1em}(2)}$$ 固定b=0.3，选择$a\in (0,1.4]$，得到Henon映射的x、y分量的分岔图如图3-4所示。

由图3-4可知，随着参数a的增加，Henon系统通过倍周期分岔进入混沌状态。对控制参数$a\in (0,1.4]$，Henon映射的Lyapunov指数图谱如图3-5所示。

由图3-5可知，当选择参数$a=1.4，b=0.3$时，Henon映射的Lyapunov指数值取到最大值且大于0，此时系统进入混沌状态，系统的相空间图和时序图分别如图3-6、3-7所示。

由图3-6和图3-7可知，当$a=1.4，b=0.3$，初始值为$\left[0,0\right]$时，Henon系统呈现经典的混沌吸引子，Henon映射的x分量和y分量的取值是随机的，系统处于混沌状态。

三、三维Lorenz连续混沌映射

Lorenz混沌系统的方程如下：$$\{\begin{array}{l}\frac{d{x}_1}{dt}=\sigma (x_2-x_1)\\ \frac{d{x}_2}{dt}=\sigma x_1-x_1{x}_3-x_2\\ \frac{d{x}_3}{dt}=x_1{x}_2-\beta x_3\end{array}\text{\hspace{1em}(3)}$$当选择参数$\sigma =10，\alpha =28，\beta =\frac{8}{3}$，初值为[1,1,1]时，系统的相空间图以及时序图分别如图3-8和图3-9所示。

由图3-8和图3-9可知，当$\sigma =10，\alpha =28，\beta =\frac{8}{3}$，初值为[1,1,1]时，Lorenz系统的三个分量的取值均是随机的，在三维空间中系统轨道围绕一点运动，显然这是一个奇怪吸引子，Lorenz系统处于混沌状态。
当$\sigma =10，\beta =\frac{8}{3}，\alpha \in (0,500]$时，Lorenz系统关于$\alpha$的分岔图如图3-10所示，最大Lyapunov指数图如图3-11所示。

由图3-11可知，当Lorenz系统的控制参数取值为$\sigma =10，\alpha =28，\beta =\frac{8}{3}$时，Lorenz系统的最大Lyapunov指数大于0，此时系统处于混沌状态。

总结

本章介绍了经典的混沌映射，包括Logistic映射、Henon映射以及Lorenz映射。通过数值仿真的方式，利用相空间图、分岔图、时序图以及Lyapunov指数图等方法，对这些经典混沌映射的动力学特征以及重要参数进行了研究。

代码

本文所有代码使用的软件均为matlab。
首先介绍计算混沌映射Lyapunov指数的方法，可以参考以下的两个网站：
链接1 Lyapunov指数计算
链接2 matlab计算Lorenz系统最大Lyapunov指数
对于一维的Logistic映射，因为只有一个Lyapunov指数，计算相对简单，可以直接采用上一篇文章中提到的一维映射Lyapunov指数计算公式进行计算（基于混沌系统的文本加密算法研究（一）——混沌及混沌加密的基础知识）。对于高维的混沌映射，由于具有两个或两个以上的Lyapunov指数，计算起来相对复杂。一般有两种方法：第一，采用上一篇文章中介绍的高维混沌映射Lyapunov指数的计算公式，即使用雅可比矩阵进行计算；第二，采用其他的数值计算方法，本文采用了链接2中提到的 G. Benettin 计算方法。
对于二维Henon映射，其具有两个Lyapunov指数，我们可以采用第一种方法，同时求解两个Lyapunov指数；对于三维Lorenz映射，因为维数较高，雅可比行列式的多次计算会导致计算结果的偏差，因此采用G.Benettin的方法求Lorenz系统的最大Lyapunov指数。
下面首先对G.Benettin的计算方法进行介绍：

首先选取系统相空间中的两个初始点，它们的距离为$d_0$（两者之间的距离足够小，可以取$10^{-10}$等），如上图$t_0$和$L(t_0)$所示。经过一次迭代，两者分别演化为$t_1$和$L^{\prime }(t_1)$，两者之间的距离发生改变。在$t_1$和$L^{\prime }(t_1)$之间的连线上，选择距离$t_1$为$d_0$的点$L(t_1)$，将$t_1$和$L(t_1)$作为新的迭代初始点。对$t_1$和$L(t_1)$进行新的一轮迭代，得到下一轮的演化结果。重复以上的连线、选择、迭代过程。直到迭代次数达到设定的次数要求。
（以上是该方法的核心思想，至于怎么计算Lyapunov指数可以参见代码部分）

下面对基于雅可比行列式的计算方法（Henon映射）进行介绍：
对给定的初始值$x_1$、$y_1$，记Henon映射的雅可比行列式为$J_1$，记第m次迭代得到的雅可比行列式为$J_m$，对乘积$J_m{J}_{m-1}\dots J_1$做QR分解：$$\begin{array}{l}qr[J_m{J}_{m-1}\cdots J_1]\\ =qr[J_m{J}_{m-1}\cdots J_2(J_1{Q}_0)]\\ =qr[J_m{J}_{m-1}\cdots J_3(J_2{Q}_1)][R_1]\\ =qr[J_m{J}_{m-1}\cdots J_4(J_3{Q}_2)][R_2{R}_1]\\ =qr[J_m{J}_{m-1}\cdots (J_i{Q}_{i-1})][R_{i-1}{R}_{i-2}\cdots R_2{R}_1]\\ =Q_m[R_m{R}_{m-1}\cdots R_2{R}_1]\\ =Q_{m}R\end{array}$$其中$qr[\bullet ]$表示QR分解。从$J_1$开始，在第i步中，我们给$Q_{i-1}$左乘矩阵$J_i$得到$B_i=J_i{Q}_{i-1}$，然后对其进行QR分解$B_i=Q_i{R}_i，i=1,2,\cdots ,m$。矩阵R是矩阵乘积$R_m\cdots R_2{R}_1$，$Q_0$为单位矩阵。由QR分解的矩阵R的性质可知，矩阵R对角线上的每个元素就是所有矩阵$R_i$的对角线上相应元素的乘积。因此，该映射的所有Lyapunov指数的近似值可以通过下式得到：$$x_k=\frac{1}{m}\sum _{i=1}^{m}ln|R_i(k,k)|,k=1,2,\cdots ,n$$对于Henon映射有两个Lyapunov指数，所以$n=2$。

1、Logistic映射

（1）分岔图

% Logistic映射的分岔图
clear
clc
mu = 0:0.001:4;                   %控制参数的取值范围与步长
x = 0.1*ones(1,length(mu));       %自变量的初始值
N1 = 1000;                        %先迭代N1次，充分迭代，排除初始值的干扰
N2 = 20;                          %将最后一次的函数值作为初始值继续进行迭代N2次并将结果作图
f = zeros(N1+N2,length(mu));      %存储迭代的函数值
for i = 1:N1+N2
    x = mu.*x.*(1-x);             %Logistic映射
    f(i,:) = x;
end
g = f(N1+1:end,:                  %使用后面的N2次迭代值作图
plot(mu,g,'k.')
xlabel('\mu')
ylabel('x')

（2）时序图

% Logistic映射的时序图
clear
clc
mu = [1,3.25,3.5,3.6];        %控制参数的取值
x = 0.4*ones(1,length(mu));   %自变量初始值
N = 500;
f = zeros(N,length(mu));      %列向量为mu取不同值时的xn迭代值
num = 1:N;                    %迭代次数向量
for i = 1:N
    x = mu.*x.*(1-x);         %Logistic映射
    f(i,:) = x;
end
for i = 1:length(mu)
    figure            
    plot(num,f(:,i),'k')      %控制参数取不同值时的logistic映射时序图
    xlabel('N')
    ylabel('x')
end

（3）Lyapunov指数——一维映射计算公式法

% Logistic映射的Lyapunov特征指数——公式
clear
clc
mu = 0:0.001:4;    %控制参数的取值
x = 0.1;           %自变量初始取值
N = 20000;
sum = log(abs(mu.*(1-2*x)));   %导数g'(x) = mu*(1-2x)
for i = 1:N-1
    x =  mu.*x.*(1-x);
    sum = sum + log(abs(mu.*(1-2*x))); %每次迭代后，进行求和
end
LE = sum/N;        %计算Lyapunov指数
plot(mu,LE,'k');
hold on
plot(mu,zeros(1,length(mu)),'k-.')     %绘制LE=0参考线
xlabel('\mu');
ylabel('Lyapunov指数')

（4）Lyapunov特征指数——G.Benettin方法

%% 最大Lyapunov指数图——G.Benettin方法
clear
clc
mu0 = 0:10^(-3):4; x0 = 0.1;
t = 200; N = 10000;
d0 = 1e-10;               
Z = [];
for j = 1:length(mu0)
    Lya = 0;
    x1 = zeros(N+t,1);          %记录两个初始点的演化轨迹
    x2 = zeros(N+t,1);
    x1(1) = x0; x2(1) = x0+d0;             %两个初始点的距离为d0
    for i = 1:N+t-1
        x1(i+1) = mu0(j)*x1(i)*(1-x1(i));     %Logistic映射，迭代一次，得到新的两个点x1和x2
        x2(i+1) = mu0(j)*x2(i)*(1-x2(i));
        d1 = abs(x1(i+1) - x2(i+1));        %计算新得到的两个点之间的距离d1
        if x2(i+1) > x1(i+1)                    %在新得到的两个点的连线之间选择距离x1为d0的点x2，继续迭代
            x2(i+1) = x1(i+1) + d0;
        else
            x2(i+1) = x1(i+1) - d0;
        end
        if i+1 > t
            Lya = Lya + log(d1/d0);     %先迭代t次，去掉过渡态，选择后面的N次来进行Lyapunov指数的计算
        end
    end
    Z = [Z; Lya/(i+1-t)];   %存储每个控制参数mu的Lyapunov指数计算值
end
plot(mu0,Z,'k')
hold on
xlabel('\mu') ylabel('Lyapunov指数')
plot(mu0,zeros(1,length(mu0)),'r-.')

可以发现，通过两种不同计算方法得到的Logistic映射Lyapunov指数图是一致的。

2、Henon映射

（1）分岔图

%% 分岔图
clear clc
a = 0:0.001:1.4;       %参数a取值范围及步长
b = 0.3;
N1 = 5000; N2 = 100;   %先迭代N1次，充分迭代后使用最后一次迭代值作为后面迭代的初始值
x = ones(N1+N2,length(a)); y = ones(N1+N2,length(a));
for j = 1:N1+N2-1
    x(j+1,:) = 1+y(j,:)-a.*x(j,:).^2;
    y(j+1,:) = b*x(j,:);
end
f = x(N1+1:end,:); g = y(N1+1:end,:);      %只用最后N2次迭代值作图
figure
plot(a,f,'k')
xlabel('a')
ylabel('x')
figure
plot(a,g,'k')
xlabel('a') 
ylabel('y')

（2）相空间图

%% 相空间图
clear clc
a = 1.4; b = 0.3;
N = 50000;
x = zeros(1,N); y = zeros(1,N);
for i = 1:N-1
    x(i+1) = 1+y(i)-a*x(i)^2;
    y(i+1) = b*x(i);
end
plot(x,y,'k.') 
xlabel('x') 
ylabel('y')

（3）时序图

%% 时序图
clear clc
M = 1000;
a = 1.4; b = 0.3;
x = zeros(1,M+1); y = zeros(1,M+1);
for i = 1:M
    x(i+1) = 1+y(i)-a*x(i)^2;
    y(i+1) = b*x(i);
end
figure
plot(0:1:M, x, 'k-') 
xlabel('n')
ylabel('x')
figure
plot(0:1:M, y, 'k-') 
xlabel('n')
ylabel('y')

（4）Lyapunov指数——雅可比行列式方法

%% Lyapunov指数图——雅可比行列式
clear clc
M = 2000;
a = 0:0.001:1.4; b = 0.3;
x = ones(M,length(a)); y = ones(M,length(a));
LE = zeros(2,length(a));           %存储参数a取不同值时对应的两个Lyapunov指数
for k = 1:length(a)
    J0 = [-2*a(k)*x(1,k), 1; b, 0];  %初始值为(a(k),b)时对应的雅可比矩阵
    Q0 = eye(2,2);                   %单位矩阵
    [Q,R] = qr(J0*Q0);               %QR分解
    Rdiag(:,1) = diag(R);            %存储矩阵R的对角线元素
    sum1 = log(abs(Rdiag(1,1)));     %求和计算Lyapunov指数
    sum2 = log(abs(Rdiag(2,1)));
    for i = 1:M-1         %迭代
        x(i+1,k) = 1 + y(i,k) - a(k)*x(i,k)^2;  %Henon映射
        y(i+1,k) = b*x(i,k);
        J = [-2*a(k)*x(i+1,k), 1; b, 0];        %迭代后的雅可比行列式
        [Q,R] = qr(J*Q);   %QR分解
        Rdiag(:,i+1) = diag(R);
        sum1 = sum1 + log(abs(Rdiag(1,i+1)));
        sum2 = sum2 + log(abs(Rdiag(2,i+1)));
    end
    LE(1,k) = sum1/M; LE(2,k) = sum2/M;          %Lyapunov指数
end
plot(a,LE(1,:),'k-.',a,LE(2,:),'k-.')
hold on
plot(a,0*ones(1,length(a)))
xlabel('a')
ylabel('Lyapunov指数')

（5）最大Lyapunov指数——G.Benettin方法

%% 最大Lyapunov指数图
clear clc
a = 0:0.001:1.4; b = 0.3;
N = 2000; t = 500;
x1 = zeros(N+t,1); y1 = zeros(N+t,1);
x2 = zeros(N+t,1); y2 = zeros(N+t,1);
d0 = 1e-10;
Z = [];
for i = 1:length(a)
    Lya = 0;
    x1(1) = 1;  y1(1) = 1;                   %初始点(x1,y1)
    x2(1) = x1(1); y2(1) = y1(1) + d0;       %初始点(x2,y2)与(x1,y1)距离为d0
    for j = 1:N+t-1
        x1(j+1) = 1 + y1(j) - a(i)*x1(j)^2; y1(j+1) = b*x1(j);
        x2(j+1) = 1 + y2(j) - a(i)*x2(j)^2; y2(j+1) = b*x2(j);
        d1 = sqrt((x1(j+1)-x2(j+1))^2+(y1(j+1)-y2(j+1))^2);    %迭代得到的新点之间的距离
        x2(j+1) = (d0/d1)*(x2(j+1)-x1(j+1))+x1(j+1);           %连线，选择距离(x1,y1)为d0的点(x2,y2)，继续迭代
        y2(j+1) = (d0/d1)*(y2(j+1)-y1(j+1))+y1(j+1);
        if j > t
            Lya = Lya + log2(d1/d0);
        end
    end
    Z = [Z Lya/(j-t)];
end
plot(a,Z,'k')
xlabel('a')
ylabel('最大Lyapunov指数')

通过两种方式得到的（最大）Lyapunov指数图如下图所示：

两个Lyapunov指数

最大Lyapunov指数

3、Lorenz映射

（1）求解Lorenz方程
可以直接使用matlab内置函数ode45求解Lorenz方程，也可以自己编写代码实现经典四级五阶Runge-Kuta法求解。
首先使用ode45求解。

%% ode45求解
clear clc
sigma = 10; alpha = 28; beta = 8/3;
x0 = ones(1,3);
[~,Y] = ode45(@(t,x) Lorenzfun(t,x,sigma,alpha,beta),[0,100],x0);   
plot3(Y(:,1), Y(:,2), Y(:,3))
xlabel('x1')
ylabel('x2')
zlabel('x3')

% Lorenz映射
function dx = Lorenzfun(t,x,sigma,alpha,beta)
dx = zeros(3,1);
dx(1) = sigma*(x(2) - x(1));
dx(2) = alpha*x(1) - x(1)*x(3) - x(2);
dx(3) = x(1)*x(2) - beta*x(3);
end

也可以自行编写程序求解。

%% 经典四级五阶龙格-库塔法求解Lorenz映射方程
clear clc
sigma = 10; alpha = 28; beta = 8/3;
x0 = 1; y0 = 1; z0 = 1;
h = 0.01; N = 10000;   %步长以及迭代次数
[x,y,z] = Runge_Kuta(x0,y0,z0,sigma,alpha,beta,h,N);
plot3(x,y,z)

配套的有四个函数，作为单独的m文件，需放在同一个文件夹里面才能运行。

% Runge-Kuta函数
function [x,y,z] = Runge_Kuta(x0, y0, z0, sigma, alpha, beta, h, N)
x = zeros(1,N); y = zeros(1,N); z = zeros(1,N);
x(1) = x0; y(1) = y0; z(1) = z0;
for j = 1:N-1
    K11 = h*RK_F(x(j), y(j), sigma);
    K21 = h*RK_G(x(j), y(j), z(j), alpha);
    K31 = h*RK_L(x(j), y(j), z(j), beta);
    K12 = h*RK_F(x(j)+K11/2, y(j)+K21/2, sigma);
    K22 = h*RK_G(x(j)+K11/2, y(j)+K21/2, z(j)+K31/2, alpha);
    K32 = h*RK_L(x(j)+K11/2, y(j)+K21/2, z(j)+K31/2, beta);
    K13 = h*RK_F(x(j)+K12/2, y(j)+K22/2, sigma);
    K23 = h*RK_G(x(j)+K12/2, y(j)+K22/2, z(j)+K32/2, alpha);
    K33 = h*RK_L(x(j)+K12/2, y(j)+K22/2, z(j)+K32/2, beta);
    K14 = h*RK_F(x(j)+K13, y(j)+K23, sigma);
    K24 = h*RK_G(x(j)+K13, y(j)+K23, z(j)+K33, alpha);
    K34 = h*RK_L(x(j)+K13, y(j)+K23, z(j)+K33, beta);
    x(j+1) = x(j)+ (K11 + 2*K12 + 2*K13 + K14)/6;
    y(j+1) = y(j)+ (K21 + 2*K22 + 2*K23 + K24)/6;
    z(j+1) = z(j) + (K31 + 2*K32 + 2*K33 + K34)/6;
end
end

function fvalue = RK_F(x, y, sigma)
fvalue = sigma*(y - x);
end

function gvalue = RK_G(x, y, z, alpha)
gvalue = alpha*x - x*z - y;
end

function lvalue = RK_L(x, y, z, beta)
lvalue = x*y - beta*z;
end

求解得到x、y、z后，即可画出相空间图以及时序图。此处略去。

（2）分岔图
使用庞加莱截面的方法画Lorenz映射的分岔图，原理如下图所示：

如图所示，用庞加莱截面$x_1=x_2$（即与$x_1-x_2$平面垂直且与$x_1、x_2$坐标轴均成45度角的平面）去截取Lorenz系统的运动轨迹。我们要求的就是Lorenz系统的轨迹与庞加莱截面的相交部分。为此，我们在$x_1-x_2-x_3$三维空间中，投影到$x_1-x_2$平面得到Lorenz系统的投影轨迹以及庞加莱截面的投影。那么我们要求的就是投影轨迹与直线$x_1=x_2$相交的部分。我们通过判断点$(x_{1,j-1},x_{2,j-1})$迭代后得到的点$(x_{1,j},x_{2,j})$，两个点是否处于直线$x_1=x_2$的不同两侧来判断投影轨迹与该直线是否有交点，即在一次迭代中轨迹是否穿过了庞加莱截面。如上图所示，若$x_{1,j-1}<x_{2,j-1}，x_{1,j}>x_{2,j}$，则两个点位于直线的两侧，通过线性插值的方法可以求出两个点的连线与直线$x_1=x_2$的交点。

%% 庞加莱截面法画分岔图
clear clc
sigma = 10; alpha = 0:1:500; beta = 8/3;
x0 = ones(1,3);
X = [];
for m = 1:length(alpha)
    [~,Y1] = ode45(@(t,x) Lorenzfun(t,x,sigma,alpha(m),beta),[0,1],x0); 
    x0 = Y1(end,:);           %先充分迭代消除初始值的影响，并使用最后一次迭代值最后下一次迭代的初始值
    [~,Y2] = ode45(@(t,x) Lorenzfun(t,x,sigma,alpha(m),beta),[0 50],x0);
    Y2(:,1) = Y2(:,2) - Y2(:,1);      %对迭代值的x1分量和x2分量作差，通过正负号交替以及线性插值法寻找满足x1=x2的点
    for j = 2:length(Y2)
        if Y2(j,1) < 0        % 即x2(j) < x1(j)
            if Y2(j-1,1) > 0  % 即x2(j-1) > x1(j-1)
                x = Y2(j,2) + Y2(j,1)*((Y2(j,2) - Y2(j-1,2))/(Y2(j-1,1) - Y2(j,1)));   
                 %在(x1(j-1), x2(j-1))以及(x1(j), x2(j))之间通过线性插值法寻找x1 = x2的点
                x = alpha(m) + abs(x)*i;         %将x作为虚部，写成复数形式，方便最后的分岔图作图
                X = [X, x];
            end
        else                  % 即x2(j) >= x1(j)
            if Y2(j-1,1) < 0  % 即x2(j-1) < x1(j-1)
                x = Y2(j,2) + Y2(j,1)*((Y2(j,2) - Y2(j-1,2))/(Y2(j-1,1) - Y2(j,1)));     %线性插值
                x = alpha(m) + abs(x)*i;
                X = [X, x];
            end
        end
    end
end
plot(X, 'k.')    %X为一个虚数组成的集合，虚数的实部为alpha的取值，虚部为对应的满足x1=x2的点的x1/x2取值
xlabel('alpha')
ylabel('x1')

% Lorenz映射
function dx = Lorenzfun(t,x,sigma,alpha,beta)
dx = zeros(3,1);
dx(1) = sigma*(x(2) - x(1));
dx(2) = alpha*x(1) - x(1)*x(3) - x(2);
dx(3) = x(1)*x(2) - beta*x(3);
end

（3）最大Lyapunov指数——G.Benettin方法

% 最大Lyapunov指数
clear clc
sigma = 10; alpha = 0:0.01:500; beta = 8/3;
d0 = 1e-7;
t = 50; n = 100;
Z = [];
for i = 1:length(alpha)
    Lya = 0;
    x0 = 1; y0 = 1; z0 = 1;
    x1 = 1; y1 = 1; z1 = z0+d0;
    for j = 1:t+n
        [T1,Y1] = ode45(@(t,x) Lorenzfun(t,x,sigma,alpha(i),beta),[0,1],[x0;y0;z0]);
        [T2,Y2] = ode45(@(t,x) Lorenzfun(t,x,sigma,alpha(i),beta),[0,1],[x1;y1;z1]);
        x0 = Y1(end,1); y0 = Y1(end,2); z0 = Y1(end,3);
        x1 = Y2(end,1); y1 = Y2(end,2); z1 = Y2(end,3);
        d1 = sqrt((x0-x1)^2+(y0-y1)^2+(z0-z1)^2);
        x1 = x0 + (d0/d1)*(x1-x0);
        y1 = y0 + (d0/d1)*(y1-y0);
        z1 = z0 + (d0/d1)*(z1-z0);
        if j > t
            Lya = Lya + log(d1/d0);
        end
    end
    Z = [Z Lya/(j-t)];
end
plot(alpha,Z,'k')
xlabel('\alpha')
ylabel('Lyapunov指数')

%% Lorenz映射
function dx = Lorenzfun(t,x,sigma,alpha,beta)
dx = zeros(3,1);
dx(1) = sigma*(x(2) - x(1));
dx(2) = alpha*x(1) - x(1)*x(3) - x(2);
dx(3) = x(1)*x(2) - beta*x(3);
end

随着维数的增加，计算Lyapunov指数所需的计算量和时间也会增加，特别是当选择较小的步长，比如控制参数a的步长较小、取值范围较大时，所需的时间更长。因此，建议先通过较大的步长计算最大Lyapunov指数，随后选择感兴趣的范围，取较小的步长进一步研究。