(完整版)基于matlab的Lorenz系统仿真研究

基于Matlab的Lorenz系统仿真研究
摘要：本文利用matlab这一数学工具对Lorenz系统进行了研究。

首先使用matlab 分析求解Lorenz方程，利用matlab的绘图功能，直观地观察了Lorenz 混沌吸引子的三维图形，并简单观察了Lorenz混沌系统对初值的敏感性；
然后对Lorenz系统进行仿真，比较分析在不同参数下的Lorenz系统仿真结果；最后验证了通过添加反馈控制的方式，可以使Lorenz方程不稳定的平衡点成为稳定的平衡点。

关键词：Lorenz系统；matlab；混沌系统
1.引言
Lorenz方程是由美国著名的气象学家Lorenz在1963年为研究气候变化，通过对对流实验的研究，建立的三个确定性一阶非线性微分方程。

这三个方程是混沌领域的经典方程，Lorenz系统也是第一个表现奇怪吸引子的连续动力系统，具有着举足轻重的作用。

Lorenz方程的表达式如下：
{
dx
dt
=σ(y−x) dy
dt
=(μ−z)x−y dz
dt
=−bz+xy
其中，σ、μ、b为正实常数。

本文利用matlab这一数学工具，对Lorenz系统进行了研究，得到了仿真结果，加深了对Lorenz系统的认识。

2.matlab求解Lorenz方程并绘图
首先建立m文件“Lorenz.m”来定义Lorenz方程，固定σ=10，μ=30，b=8/3，程序如下所示：
function dx=Lorenz(t,x)
dx=[-10*(x(1)-x(2));30*x(1)-x(2)-x(1)*x(3);x(1)*x(2)-2.6667*x(3)];
end
然后利用ode45命令来求解Lorenz方程并绘制图形，初值取x=y=z=0.1。

程序如下所示：
>> clf
>> x0=[0.1,0.1,0.1];
>> [t,x]=ode45('Lorenz',[0,100],x0);
>> subplot(2,2,1)
>> plot(x(:,1),x(:,3))
>> title('(a)')
>> subplot(2,2,2)
>> plot(x(:,2),x(:,3))
>> title('(b)') >> subplot(2,2,3)
>> plot(x(:,1),x(:,2)) >> title('(c)') >> subplot(2,2,4)
>> plot3(x(:,1),x(:,2),x(:,3)) >> title('(d)')
运行上述程序，可得到如下波形：
其中，图（a ）为Lorenz 混沌吸引子在x-z 平面上的投影，图（b ）为Lorenz 混沌吸引子在y-z 平面上的投影，图（c ）为Lorenz 混沌吸引子在x-y 平面上的投影，图（d ）为Lorenz 混沌吸引子的三维图。

可以看到，混沌吸引子在各平面上的投影类似于横写的“8”字形。

由于参数σ=10，μ=30，b=8/3时为混沌系统，对初值具有敏感性，初值很小的差异会引起系统行为的显著改变。

因此，将初值改为x=z=0.1,y=0.11,绘制此时混沌吸引子在x-z 平面上的投影，并与初值为x=y=z=0.1时混沌吸引子在x-z 平面上的投影放在同一张图中比较。

为了区别两者，初值为x=y=z=0.1时混沌吸引子在x-z 平面上的投影用蓝色，初值改为x=z=0.1,y=0.11时混沌吸引子在x-z 平面上的投影用红色。

程序如下所示：
>> clf
>> x0=[0.1,0.1,0.1];
>> [t,x]=ode45('Lorenz',[0,100],x0); >> plot(x(:,1),x(:,3)) >> hold on
>> x0=[0.1,0.1,0.1];
>> [t,x]=ode45('Lorenz',[0,100],x0); >> x0=[0.1,0.11,0.1];
(a)
(b)
(c)
(d)
>> [t,x]=ode45('Lorenz',[0,100],x0); >> plot(x(:,1),x(:,3),'r*')
得到的图形如下所示：
可以看到，虽然初值只有0.01的改变，红色与蓝色图形明显不重合，这证明了系统的敏感性。

3.matlab 对Lorenz 系统仿真
首先利用matlab 的Simulink 功能，搭建Lorenz 系统的模型，仿真模型如下图所示：
010
20
30
40
50
60
在仿真模型中，取参数σ=10，b=8/3，观察参数μ取不同值时系统的运行状态。

根据文献[1]的分析，当参数0<μ<1时，只有一个稳定平衡点O（0,0,0）。

取初值为x=y=z=2,参数μ=0.5，仿真停止时间取为50，运行仿真。

得到x、y、z 的相图以及x-z，y-z，x-y的图形依次如下所示：
050100150 0
0.5
1
1.5
2
2.5
3
可见，系统很快地趋向并稳定在O （0,0,0），验证了前面所述。

根据文献[1]，当μ>1时，系统有三个平衡点：原点O(0,0,0)和P+，P-。

此时原点的特征值中有正值，因此原点为鞍点，是不稳定平衡点。

当1<μ<13.926时，不稳定流形最终螺旋地趋于与之同侧的平衡点P+或P-；当μ=13.926时，不稳定流形刚好无限趋于原点O ，即出现同宿轨；当μ>13.926时，不稳定流形将绕到另一侧，最终趋于与之异侧的P+或P-。

可见，μ是一个同宿分岔点。

因此，取初值x=y=z=2，μ=8，仿真停止时间为50，运行仿真，得到x 、y 、z 的相图以及x-z ，y-z ，x-y 的图形依次如下所示：
050100150200250300350
2
3456
7150200250300350400450
2
3456
7
100200300400500600700
24681012
可以看到，系统趋于与之同侧的平衡点P+或P-。

取初值x=y=z=2，μ=18，仿真停止时间为50，运行仿真，得到x 、y 、z 的相图以及x-z ，y-z ，x-y 的图形依次如下所示：
050010001500
-20
-15-10-50510
15
050010001500
10
15
20
25
30
可以看到，系统趋于与之同侧的平衡点P+或P-。

为了观察μ=13.926的同宿分岔点现象，在μ=13.926附近不断尝试，最终在μ=15.39682328时观察到比较明显的过渡迹象。

取初值x=y=z=2，μ=15.39682328，仿真停止时间为50，运行仿真，得到x、y、z的相图以及x-z，y-z，x-y的图形
依次如下所示：
050010001500
-10
-5
5
10
15
20
500
1000
1500
-15-10
-5
5
10
15
500
1000
1500
2000
2500
3000
051015202530354045
50
可以看到，虽然最终轨线趋向于与之同侧的平衡点P+或P-，但有着明显的过渡迹象。

可以推测，当μ取15.39682328到15.39682330间的某一个数值时，会出现同宿轨现象。

根据文献[1]，当μ>24.74时，P+与P-变为不稳定的，也就是说系统进入“混沌区”。

此时三个平衡点O、P+、P-都不稳定。

取初值x=y=z=2，μ=30，仿真停止时间为100，运行仿真，得到x、y、z的相图以及x-z，y-z，x-y的图形依次如下所示：
0500100015002000250030003500400045005000
-20
20
40
60
80
100
0500100015002000250030003500400045005000
-40
-20
20
40
60
80
100
500
1000
1500
2000
2500
3000
3500
4000
4500
5000
0102030405060708090
100
可以看到，上述图形中，轨线绕着P+若干圈后，又绕着P-若干圈，如此循环，符合文献[1]的描述。

为了观察由系统趋向于与之异侧的平衡点向系统的混沌状态的过渡现象，在μ=24.74附近反复不断尝试，最终发现当μ=23.299时，可以观察到明显的过渡迹象。

因此，取初值x=y=z=2，μ=23.299，仿真停止时间为100，运行仿真，得到x、y、z的相图以及x-z，y-z，x-y的图形依次如下所示：
5001000150020002500
-20
-10
10
20
30
40
05001000150020002500
-20
-10
10
20
30
40
500
1000
1500
2000
2500
051015202530354045
50
可以看到，在上图中，轨线看起来稳定在一条围绕与之异侧的平衡点的轨道上。

仅从仿真运行的这段时间，无法判断系统是处于混沌状态还是会趋向于与之异侧的平衡点，可以看出明显的过渡迹象。

4.对Lorenz 系统的反馈控制
系统的稳定是系统的基本要求。

为了使Lorenz 系统的不稳定平衡点变为稳定平衡点，根据文献[2]，可以通过加入反馈控制的方法实现。

加入反馈后，Lorenz 方程变为：
{ dx dt =σ(y −x)dy dt =(μ−z )x −y dz dt =−bz +xy
+ky 由上式可以看出，第二个方程加入了简单的线性反馈ky 。

建立加入反馈后的系统仿真模型，如下图所示：
根据文献[2]的分析，当k<-35.1时，可以满足使系统稳定的要求。

取初值x=y=z=2，μ=30，仿真停止时间为100，k=-36，运行仿真，得到x 、y 、z 的相图以及x-z ，y-z ，x-y 的图形依次如下所示：
0100200300400500600
02
4
6
8
10
12
020040060080010001200
05
10
15
20
25
0100200300400500600
02
4
6
8
10
12
可以看到，系统很快趋于原点O(0,0,0)并稳定下来，这验证了通过加入反馈使Lorenz系统变得稳定的这一方法的正确性。

5.结论
本文直观地观察了Lorenz混沌吸引子的三维图形，并简单观察了Lorenz混
沌系统对初值的敏感性，比较分析了在不同参数下的Lorenz系统仿真结果，最后
验证了添加反馈控制这一方法可以使Lorenz方程不稳定的平衡点成为稳定的平
衡点。

通过使用matlab对Lorenz系统仿真，直观地观察到了Lorenz系统的运行
轨迹，加深了对Lorenz方程和混沌现象的理解。

参考文献：
[1]刘崇新.非线性电路理论及应用[M].西安：西安交通大学出版,2007.201-208
[2]朱少平.Lorenz方程的动力学特性与控制[J].陕西教育学院学报，2007,23（4）：
81-84
[3]赖宏慧.基于matlab的Lorenz系统模拟实验仿真[J].科技信息.2010,17:18-19
[4]柴彩春.关于Lorenz方程的动力学性态研究[J].廊坊师范学院学报，2011,11（1）：8-10
[5]刘庆花.基于Matlab的lorenz混沌系统仿真[J].现代商贸工业，2014,02:194-195。