Homework10

李尧 2013301020048 天眷班
problem3.26

混沌 洛伦兹模型 相图

摘要
洛伦兹模型是研究混沌的经典的模型之一。本文中将会以他为例分析混沌的行为，并以图片的形式给出一定的初始条件下的数值的解。此外，本文中将会尝试使用数值计算的方法处理较简单的桌球的混沌问题。
关键字：混沌 洛伦兹模型 相图 桌球

背景：
在1972年12月29日，美国麻省理工教授、混沌学开创人之一E.N.洛仑兹在美国科学发展学会第139次会议上发表了题为《蝴蝶效应》的论文，提出一个貌似荒谬的论断：在巴西一只蝴蝶翅膀的拍打能在美国得克萨斯州产生一个龙卷风，并由此提出了天气的不可确预报性。至此以后，人们对于混沌学研究的兴趣十分浓厚。其中洛伦兹提出的用于描述天气变化的非线性方程组被称之为洛伦兹模型，是一个典型的混沌系统。本文将就理论和数值解分析这一系统。

物理数学模型和数值计算方法
原始的洛伦兹的三个方程如下：
$\frac{dx}{dt}=\sigma(y-x)$
$\frac{dy}{dt}=-xz+rx-y$
$\frac{dz}{dt}=xy-bz$
其中想x,y和z分别表示温度，密度和速度，其他都是待定的系数。为了方便研究，我们取$\sigma=10$,$b=\frac{8}{3}$,x=1，y=0,z=0.只考虑r对整个方程组的解的影响。具体的计算方法采用欧拉法，对应地把上式改写成迭代的形式：
$x_{i+1}=x_i+\Delta t\sigma(y_i-x_i)$
$y_{i+1}=y_i+\Delta t(-x_iz_i+rx_i-y_i)$
$z_{i+1}=z_i+\Delta t(x_iy_i-bz_i)$
即可数值求解。
代码详见program

数值解稳定性的讨论
考虑在数值解附近有 $\delta x$,$\delta y$,$\delta z$的偏离,代入以上方程只保留一阶项得：
$\frac{d}{dt}\delta x=\sigma*(\delta y-\delta x)$
$\frac{d}{dt}\delta y=-x\delta -z\delta x-\delta y+r\delta x$
$\frac{d}{dt}\delta z=x\delta y+y\delta x-b\delta z$
改写成矩阵的形式
$\frac{d}{dt}\delta \vec{r}=M(r)\delta \vec{r}$
若M的每一个特征值的实部都小于零，那么解就稳定，否则不稳定。这里我们就本文中使用的初始条件下的解的稳定性做简单的讨论。
当r<=1时，原方程只有一个稳定点（0,0,0），代入上式得：
$\lambda_1*\lambda_2=10-10*r \lambda_1+\lambda2=-11$
也就是说必然不会出现混沌。当r>1时，原方程有三个稳定点，分别是（0,0,0），（$(b(r-1))^{\frac{1}{2}},(b(r-1))^{\frac{1}{2}},r-1$）和($-(b(r-1))^{\frac{1}{2}},-(b(r-1))^{\frac{1}{2}},r-1$).其中（0,0,0）已经讨论过了,结论是当r>=1时不稳定。可以证明不同的稳定点的稳定条件是不会重合的。带入另外两个稳定点，可得一个大于1小于$\frac{41}{3}$的范围.这和混沌的边界条件不符。
为什麽会这样呢？我们应该把上面的微分方程改成迭代的离散形式。只要迭代的$\Delta t$趋于零，那么数值算法本身的稳定性就对结果没有影响。这时，我们可以定义n阶的稳定点。如果点$r_i$代入得到$r_(i+1)$,即
$r_{i+1}=f(r_i)$
且满足 $r_i=r_{n+i}$的条件，我们就称点$r_i$为n阶的稳定点。前面的讨论中我们实际上只计算了1阶的情况。周期性的解从离散的角度来看就是n很大的稳定点。
求解方程
$r_i=f^n(r_i)$
得到所有的稳定点，并按以上的方法逐个讨论稳定条件，就可以得到混沌的条件。实际上就洛伦兹模型而言只研究n趋于无穷的情况就行了，但是对于一个非线性的方程而言，要解出n解得稳定点十分困难！
关于（0,0,0）的结论会在下文中被计算机验证。
当r=0.9时：

当r=1.0时：

当r=1.1时：

可见在临界值附近稳定的解有一个跳跃。

程序运行结果
试运行程序，绘制r=25时的t-z,x-z,x-z(y=0)和y-z(x=0):
picture1:

picture2:

picture3:

picture4:

接下来改变参数，尝试r=20，相应的相图为：
picture5：

y-z，x=0和x-z，y=0不存在。
虽然理论的临界点实如24.72，但是由于计算误差等原因，实际观测到在r=20.9附近就已经变成混沌的。
picture6：

picture7：

picture8：