产品
摘要：混沌（Chaos）是指发生在确定系统中的貌似随机的不规则运动，本文将基于经典的二维系统，然后根据动力学方程研究其混沌产生过程以及相对应的MATLAB仿真，再讨论Lyapunov指数以及正平衡点。
上一篇中介绍了一维系统，这次我们将维数提升到二。
/*仅当作学习笔记，若有纰漏欢迎友好交流指正，此外若能提供一点帮助将会十分荣幸*/
0引言
对于一维离散混沌系统来说，其具有形式简单、产生混沌序列时间短等优点。但其缺点是空间太小，难以处理较为复杂的情况。因此本文将从经典的二维混沌系统入手以此为基础，研究动力学系统的发展、混沌的产生以及相应的MATLAB仿真，借此建立一个对二维混沌系统有更为完整的认识与理解。
1. 厄农（Henon）映射
1.1系统的混沌
厄农映射是由法国数学家米歇尔·厄农提出的，以此作为洛伦茨模型的庞加莱截面的简化模型。对经典厄农映射而言，任意初始点或趋向厄农奇异吸引子，或发散至无穷大。厄农吸引子具有分形结构，其在一个方向上连续，另一个方向上则为一个康托尔集，其方程动力学微分方程为：
其中：a、b为系统的参数，分别取为0.3、0.4，其对应MATLAB程序如下：
%厄农映射
%x(j+1)=1+yj-a*xj*xj
%y(j+1)=b*xj
clc
clear all
close all
a=0.4;b=0.3;%系统参数设置
xjL=[];yjL=[];
for k=1:140%迭代次数
a=0.5+(k-1)*0.01;
xj=[];yj=[];
x=0.1;y=0.1;
for j=1:200
x=1+y-a*x*x;y=b*x;%系统函数
xj=[xj;j x];
yj=[yj;j y];
if j>=120
xjL=[xjL;a b x];
yjL=[yjL;a b y];
end
end
end
plot(xjL(:,1),xjL(:,3),'k.')
%坐标轴设置
xlabel('\itt','FontSize',20,'FontName','Times New Roman');
ylabel('\itx','FontSize',20,'FontName','Times New Roman');
grid on;仿真结果：
其中，横轴代表迭代时间，纵轴为迭代值x.从图中可以看出，在出现两次分岔之后，系统进入混沌状态混沌状况。
1.2 最大Lyapunov指数
根据系统参数a的变化，计算得到函数的最大Lyapunov指数，其程序如下：
clc, clear
a=0:0.01:1.6; n=140; S=[];
for j=1:length(a)
x=0.2; y=0.3; %x,y的初值
for i=1:n
x2=1-a(j)*x^2+y; y=0.3*x; x=x2; %首先进行序列迭代
end
if x(end)>-100 & x(end)<100 %若不发散，再计算指数
x=0.2; y=0.3; %设立二阶单位矩阵
JJ=eye(2);
for i=1:n
x2=1-a(j)*x^2+y; y=0.3*x; x=x2;
J=[-2*a(j)*x, 1; 0.3, 0]; %分别对两迭代式对于x、y求偏导数
JJ=JJ*J;
end
L=eigs(JJ,1); %求模最大的特征值
S=[S,[a(j);log(abs(L))/n]]; %把a值及指数值保存
end
end
plot(S(1,:),S(2,:),'.-') %画出a对应的最大Lyapunov指数
hold on, plot([a(1),a(end)],[0,0],'k')
xlabel('\it a'), ylabel('最大Lyapunov指数')得最大Lyapunov指数变化图：
1.3系统正平衡点的求解
根据微分方程（1），利用平衡条件：
可得：
变化一下，可得：
利用matlab的solve函数求解，其程序为：
syms S L
a=0.4;b=0.3;%系统参数设置
[solS,solL]=solve(1-a*S^2+L-S==0,b*S-L==0,S,L)%列出方程组
solutions=[solS,solL]
x1=double(solS)%an*的值
x2=double(solL)%bn*的值计算结果为：
可得正平衡点E*为：
（0.9321，0.2796）
2. Lotka-Volterra捕食模型的混沌
在捕食系统中，最经典的捕食模型就是Lotka-Volterra模型(简称L-V模型)，其模型方程为：
L-V模型中，、
分别表示在t时刻食饵和捕食者的密度,
、
则分别表示食饵和捕食者在t时刻的种群密度变化量,a10代表食饵种群内自然增长率、a12代表捕食者对食饵捕食的所造成的影响（或者说威胁）、a20代表捕食者种群内部因素带来的数量变化，比如生老病死亦或同族之间的相互争斗带来的影响、a21代表食饵对捕食者所提供的积极影响，即提供食物来源。
2.1系统的混沌发展
把系统参数中的a10当作一个变化的参数，而另外的a12、a20、a21当作常数，就可以利用类似于得到Logitic模型分叉图那样调制得到L-V模型的系统分叉图，其相应程序为：
clc
clear all
close all
a10=1.4;
xjL=[];yjL=[];
for k=1:120
a10=1.4+(k-1)*0.001;
xj=[];yj=[];
x=0.6667,y=3.500;
for j=1:20
x=x*(a10-0.4*y);
y=y*(-0.2+0.3*x);
xj=[xj;j x];
yj=[yj;j y];
if j>=15
xjL=[xjL;a10 x];
yjL=[yjL;a10 y];
end
end
end
plot(xjL(:,1),xjL(:,2),'.k')
xlabel('\itt','FontSize',20,'FontName','Times New Roman');
ylabel('\itx','FontSize',20,'FontName','Times New Roman');
grid on;得分岔图：
从上图可以看到，从起点出发，随着迭代次数的增加（即繁殖代数），物种x1开始慢慢出现出现6条分岔。
2.2最大Lyapunov指数
根据系统参数a10的变化，计算得到函数的最大Lyapunov指数，其程序如下：
clc, clear
a=0:0.01:1.6; n=140; S=[];
for j=1:length(a)
x=0.6667,y=3.500; %x,y的初值
for i=1:n
x2=a(j)*x-0.4*x*y; y2=1.2*y+0.3*x*y;x=x2; %首先进行序列迭代
end
if x(end)>-100 & x(end)<100 %若不发散，再计算指数
x=0.6667,y=3.500;
JJ=eye(2);
for i=1:n
x2=a(j)*x-0.4*x*y; y2=1.2*y+0.3*x*y; x=x2;
J=[a(j)-0.4*y,-0.4*x;0.3*y, 1.2+0.3*x];
JJ=JJ*J;
end
L=eigs(JJ,1); %求模最大的特征值
S=[S,[a(j);log(abs(L))/n]]; %把a值及指数值保存
end
end
plot(S(1,:),S(2,:),'.-') %画出a对应的最大Lyapunov指数
hold on, plot([a(1),a(end)],[0,0],'k')
xlabel('\it a'), ylabel('最大Lyapunov指数')仿真结果为：
2.3正平衡点的求解
我们设平衡点E*（x1*,x2*）为由于各个种群的变化率为0，由上文中的模型可得；
求解该方程组，我们可以用到solve函数来求解，其对应程序如下：
%计算正平衡点
symsxy%定义自变量
symsa10a12a20a22a21%定义系数
[solx,soly]=solve(x*(a10-a12*y)==0,y*(-a20+a21*x)==0,x,y)
solutions=[solx,soly]%方程组
m=double(solx)%x1*的值
k=double(soly)%x2*的值结果为：
带入系统参数：
可求得正平衡点E*（0.6667，3.500）
3. 总结
根据上述例子，我们可以看到对于二维系统而言，当有两个系统参数时，为了简化计算我们可以尝试着只改变其中一个参数，而将另外一个参数设为常数。当然对于多个系统参数的系统，也可以只设置一个变量参数，另外几个设置为常量参数。当然如果需要讨论不同参数的作用时，我们可以采用控制变量法，依次讨论各个参数的迭代变化所带来的影响。
最大Lyapunov指数的求取也可以对系统的混沌状态进行刻画，即大于0时，系统是处于混沌状态。而正平衡点的求取，可以找出系统所处平衡状态的节点，这也对找出系统的混沌，以及对混沌状态的模拟提供了依据。
参考文献
[1] 柯赟,蒋天发.基于离散余弦和Contourlet混合变换域的图像水印方案.武汉大学学报：工学版,2012,45(2):264-267.
155-2731-8020
