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

  • 前言
  • 一、一维Logistic混沌映射
  • 二、二维Henon混沌映射
  • 三、三维Lorenz连续混沌映射
  • 总结
  • 代码
    • 1、Logistic映射
    • 2、Henon映射
    • 3、Lorenz映射

前言

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

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

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

一、一维Logistic混沌映射

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

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

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

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

二、二维Henon混沌映射

Henon映射是一个经典的二维离散混沌映射。其方程如下： { x n + 1 = 1 + y n − a x n 2 y n + 1 = b x n (2) \left\{ \begin{array}{l} x_{n+1}=1+y_n-ax_n^2\\ y_{n+1}=bx_n \end{array}\right. \tag{2} {xn+1​=1+yn​−axn2​yn+1​=bxn​​(2) 固定b=0.3，选择 a ∈ ( 0 , 1.4 ] a\in(0,1.4] a∈(0,1.4]，得到Henon映射的x、y分量的分岔图如图3-4所示。

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

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

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

三、三维Lorenz连续混沌映射

Lorenz混沌系统的方程如下： { d x 1 d t = σ ( x 2 − x 1 ) d x 2 d t = σ x 1 − x 1 x 3 − x 2 d x 3 d t = x 1 x 2 − β x 3 (3) \left\{ \begin{array}{l} \frac{dx_1}{dt}=\sigma(x_2-x_1)\\ \frac{dx_2}{dt}=\sigma x_1-x_1x_3-x_2\\ \frac{dx_3}{dt}=x_1x_2-\beta x_3 \end{array}\right. \tag{3} ⎩ ⎨ ⎧​dtdx1​​=σ(x2​−x1​)dtdx2​​=σx1​−x1​x3​−x2​dtdx3​​=x1​x2​−βx3​​(3)当选择参数 σ = 10 ， α = 28 ， β = 8 3 \sigma=10，\alpha=28，\beta=\frac{8}{3} σ=10，α=28，β=38​，初值为[1,1,1]时，系统的相空间图以及时序图分别如图3-8和图3-9所示。


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

由图3-11可知，当Lorenz系统的控制参数取值为 σ = 10 ， α = 28 ， β = 8 3 \sigma=10，\alpha=28，\beta=\frac{8}{3} σ=10，α=28，β=38​时，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 d_{0} d0​（两者之间的距离足够小，可以取 1 0 − 10 10^{-10} 10−10等），如上图 t 0 t_{0} t0​和 L ( t 0 ) L(t_{0}) L(t0​)所示。经过一次迭代，两者分别演化为 t 1 t_{1} t1​和 L ′ ( t 1 ) L'(t_{1}) L′(t1​)，两者之间的距离发生改变。在 t 1 t_{1} t1​和 L ′ ( t 1 ) L'(t_{1}) L′(t1​)之间的连线上，选择距离 t 1 t_{1} t1​为 d 0 d_{0} d0​的点 L ( t 1 ) L(t_{1}) L(t1​)，将 t 1 t_{1} t1​和 L ( t 1 ) L(t_{1}) L(t1​)作为新的迭代初始点。对 t 1 t_{1} t1​和 L ( t 1 ) L(t_{1}) L(t1​)进行新的一轮迭代，得到下一轮的演化结果。重复以上的连线、选择、迭代过程。直到迭代次数达到设定的次数要求。
（以上是该方法的核心思想，至于怎么计算Lyapunov指数可以参见代码部分）

下面对基于雅可比行列式的计算方法（Henon映射）进行介绍：
对给定的初始值 x 1 x_{1} x1​、 y 1 y_{1} y1​，记Henon映射的雅可比行列式为 J 1 J_{1} J1​，记第m次迭代得到的雅可比行列式为 J m J_{m} Jm​，对乘积 J m J m − 1 … J 1 J_{m}J_{m-1}\ldots J_{1} Jm​Jm−1​…J1​做QR分解： q r [ J m J m − 1 ⋯ J 1 ] = q r [ J m J m − 1 ⋯ J 2 ( J 1 Q 0 ) ] = q r [ J m J m − 1 ⋯ J 3 ( J 2 Q 1 ) ] [ R 1 ] = q r [ J m J m − 1 ⋯ J 4 ( J 3 Q 2 ) ] [ R 2 R 1 ] = q r [ J m J m − 1 ⋯ ( J i Q i − 1 ) ] [ R i − 1 R i − 2 ⋯ R 2 R 1 ] = Q m [ R m R m − 1 ⋯ R 2 R 1 ] = Q m R \begin{array}{l} qr[J_mJ_{m-1}\cdots J_1]\\=qr[J_mJ_{m-1}\cdots J_2(J_1Q_0)]\\ =qr[J_mJ_{m-1}\cdots J_3(J_2Q_1)][R_1]\\ =qr[J_mJ_{m-1}\cdots J_4(J_3Q_2)][R_2R_1]\\=qr[J_mJ_{m-1}\cdots (J_iQ_{i-1})][R_{i-1}R_{i-2}\cdots R_2R_1]\\=Q_m[R_mR_{m-1}\cdots R_2R_1]\\=Q_mR\end{array} qr[Jm​Jm−1​⋯J1​]=qr[Jm​Jm−1​⋯J2​(J1​Q0​)]=qr[Jm​Jm−1​⋯J3​(J2​Q1​)][R1​]=qr[Jm​Jm−1​⋯J4​(J3​Q2​)][R2​R1​]=qr[Jm​Jm−1​⋯(Ji​Qi−1​)][Ri−1​Ri−2​⋯R2​R1​]=Qm​[Rm​Rm−1​⋯R2​R1​]=Qm​R​其中 q r [ ∙ ] qr[\bullet] qr[∙]表示QR分解。从 J 1 J_1 J1​开始，在第i步中，我们给 Q i − 1 Q_{i-1} Qi−1​左乘矩阵 J i J_i Ji​得到 B i = J i Q i − 1 B_i=J_iQ_{i-1} Bi​=Ji​Qi−1​，然后对其进行QR分解 B i = Q i R i ， i = 1 , 2 , ⋯ , m B_i=Q_iR_i，i=1,2,\cdots,m Bi​=Qi​Ri​，i=1,2,⋯,m。矩阵R是矩阵乘积 R m ⋯ R 2 R 1 R_m\cdots R_2R_1 Rm​⋯R2​R1​， Q 0 Q_0 Q0​为单位矩阵。由QR分解的矩阵R的性质可知，矩阵R对角线上的每个元素就是所有矩阵 R i R_i Ri​的对角线上相应元素的乘积。因此，该映射的所有Lyapunov指数的近似值可以通过下式得到： x k = 1 m ∑ i = 1 m l n ∣ R i ( k , k ) ∣ , k = 1 , 2 , ⋯ , n x_k=\frac{1}{m}\sum_{i=1}^mln|R_i(k,k)|,k=1,2,\cdots,n xk​=m1​i=1∑m​ln∣Ri​(k,k)∣,k=1,2,⋯,n对于Henon映射有两个Lyapunov指数，所以 n = 2 n=2 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+N2x = 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:Nx = 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-1x =  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;             %两个初始点的距离为d0for i = 1:N+t-1x1(i+1) = mu0(j)*x1(i)*(1-x1(i));     %Logistic映射，迭代一次，得到新的两个点x1和x2x2(i+1) = mu0(j)*x2(i)*(1-x2(i));d1 = abs(x1(i+1) - x2(i+1));        %计算新得到的两个点之间的距离d1if x2(i+1) > x1(i+1)                    %在新得到的两个点的连线之间选择距离x1为d0的点x2，继续迭代x2(i+1) = x1(i+1) + d0;elsex2(i+1) = x1(i+1) - d0;endif i+1 > tLya = Lya + log(d1/d0);     %先迭代t次，去掉过渡态，选择后面的N次来进行Lyapunov指数的计算endendZ = [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-1x(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-1x(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:Mx(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)));endLE(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)距离为d0for j = 1:N+t-1x1(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 > tLya = Lya + log2(d1/d0);endendZ = [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-1K11 = 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 x1​=x2​（即与 x 1 − x 2 x_1-x_2 x1​−x2​平面垂直且与 x 1 、 x 2 x_1、x_2 x1​、x2​坐标轴均成45度角的平面）去截取Lorenz系统的运动轨迹。我们要求的就是Lorenz系统的轨迹与庞加莱截面的相交部分。为此，我们在 x 1 − x 2 − x 3 x_1-x_2-x_3 x1​−x2​−x3​三维空间中，投影到 x 1 − x 2 x_1-x_2 x1​−x2​平面得到Lorenz系统的投影轨迹以及庞加莱截面的投影。那么我们要求的就是投影轨迹与直线 x 1 = x 2 x_1=x_2 x1​=x2​相交的部分。我们通过判断点 ( x 1 , j − 1 , x 2 , j − 1 ) (x_{1,j-1},x_{2,j-1}) (x1,j−1​,x2,j−1​)迭代后得到的点 ( x 1 , j , x 2 , j ) (x_{1,j},x_{2,j}) (x1,j​,x2,j​)，两个点是否处于直线 x 1 = x 2 x_1=x_2 x1​=x2​的不同两侧来判断投影轨迹与该直线是否有交点，即在一次迭代中轨迹是否穿过了庞加莱截面。如上图所示，若 x 1 , j − 1 < x 2 , j − 1 ， x 1 , j > x 2 , j x_{1,j-1}<x_{2,j-1}，x_{1,j}>x_{2,j} x1,j−1​<x2,j−1​，x1,j​>x2,j​，则两个点位于直线的两侧，通过线性插值的方法可以求出两个点的连线与直线 x 1 = x 2 x_1=x_2 x1​=x2​的交点。


%% 庞加莱截面法画分岔图
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];endelse                  % 即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];endendend
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 > tLya = Lya + log(d1/d0);endendZ = [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指数，随后选择感兴趣的范围，取较小的步长进一步研究。

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

  1. 基于混沌系统的文本加密算法研究系列

    基于混沌系统的文本加密算法研究(三) 前言 一. Hodgkin-Huxley模型的数学形式 2.Hodgkin-Huxley模型的混沌分析 (1)外部电流为混沌序列 (2)外部电流为周期性刺激电流 ...

  2. 基于混沌映射的文本加密算法研究系列

    基于混沌映射的文本加密算法研究(四) 前言 一.传统DES密码算法 二.典型的文本混沌加密算法 1.Logistic映射 2.Henon映射 3.Lorenz映射 4.Hodgkin-Huxley模型 ...

  3. 混沌图像加密matlab,基于复合混沌系统的彩色图像加密算法及Matlab实现

    第27卷 第3期 湖 南 城 市 学 院 学 报 (自然科学版) Vol. 27 No.3 2018年5月 Journal of Hunan City University (Natural Scie ...

  4. 阿诺德图像加密c语言,基于Arnold置乱的数字图像加密算法（二）

    前文我们介绍了基于Arnold置乱的数字图像加密算法的两种图像置乱变换,今天我们介绍的是另外三种图像置乱变换:基于骑士巡游的图像置乱变换.基于Arnold变换的数字图像置乱和基于仿射变换的置乱变换. ...

  5. r语言 svm 大样本_r语言基于SVM模型的文本分类研究 附数据代码

    1 Perceptron 与 SVM 概念介绍 1.1 感知机 (Perceptron) 感知机( perceptron ) 1957 年由 Rosenblatt 提出,是神经网络与支持向 量机的基础 ...

  6. 基于bert模型的文本分类研究：“Predict the Happiness”挑战

    1. 前言 在2018年10月,Google发布了新的语言表示模型BERT-"Bidirectional Encoder Representations from Transformers& ...

  7. 基于BERT模型的文本分类研究 TensorFlow2实现（内附源码）【自然语言处理NLP-100例】

  8. 基于Duffing系统的分数阶混沌研究【基于matlab的动力学模型学习笔记_5】

    /*仅当作学习笔记,若有纰漏欢迎友好交流指正,此外若能提供一点帮助将会十分荣幸*/ 前面的几篇博文我们提到提到的都是整数阶模型,这里我们将对分数阶模型进行一个简单的研究. 摘要:与整数阶混沌相比,分数 ...

  9. 【图像加密】基于混沌系统进行灰度图像加密附Matlab代码

    1 简介 ​ 1 基于混沌系统的图像加密解密 Logistic混沌置乱,先不说有多复杂,其实很简单. Logistic函数是源于一个人口统计的动力学系统,其系统方程形式如下: **X(k+1) = u ...

最新文章

  1. 2022-2028年中国分散式风电行业投资分析及前景预测报告
  2. sys.check_constraints
  3. Ubuntu Vim YouCompleteMe 安装
  4. 【javaweb】Session原理以及浏览器禁止Cookie之后服务器如何获取Session
  5. File.separator或File.pathSeparator
  6. 自制Android相机
  7. k8s设置标签禁止istio边车sidebar注入
  8. TCP的困境与解决方案
  9. ubuntu 远程桌面
  10. C语言输入函数换行符赋给变量B,C语言程序设计第3章顺序结构程序设计.pptx-资源下载在线文库www.lddoc.cn...
  11. Hadoop HIVE 条件控制函数
  12. 文档比较比对工具Beyond Compare
  13. parceljs 中文文档24小时诞生记
  14. ITIL4培训系列之变更支持流程和实践讲解
  15. zemax设置 像方远心_ZEMAX|如何翻转整个光学系统
  16. 无线局域服务器架设方法,技巧：如何实现局域网架设BT服务器
  17. 笔记本处理器排名_上半年最受欢迎处理器TOP10榜单：AMD终进榜，9代酷睿无缘前10...
  18. M40Z-025003TB0西克光电开关 订货号: 1200128
  19. 社区发现(社团检测)模块度Modularity详细介绍
  20. python好看图案的编程代码_利用Python绘制了一些有意思的图案

热门文章

  1. 2--STM32+USB移植+HID 与AUDIO类MIDI设备组成的复合设备
  2. 戴尔便携式计算机右键,戴尔最新笔记本如何在右键菜单添加“在此处打开命令窗口”设置项?...
  3. 基于Simulink的步进电机仿真实现（文末资源）
  4. 计算机视觉最全专栏教程总结
  5. 关于titanic数据集（一）
  6. selenium下对指定元素进行截图
  7. 如何修改SVN的地址
  8. 网上赚钱第一步是借鉴与学习
  9. ubuntu 安装小企鹅拼音输入法
  10. 1.01.21盒子模型，浮动，定位