10 第十章 关于混沌
柴门闻犬吠，风雪夜归人。
10.1 不确定的世界：混沌
世界是可以被预测的吗？或者换句话说，我们生活的空间，这片宇宙能否用一组数学方程来描述？无论这些方程多么地复杂，或者是拥有比全宇宙分子原子总数都还要多的变量……
如果有一天我们把宇宙的方程写了出来，并且知道了宇宙诞生那一刻所有的初始条件，那是不是过去、现在和未来所有已经发生过或即将要发生的事情，我们都能计算出来？我们是否生活在一个可计算的宇宙里？这个问题被称为“拉普拉斯妖”，拉普拉斯妖是存在于哲学尽头的一只妖怪，如果它是真实的，那人类所谓的自由意志是不是从一开始就被固定住了呢？
于我自己而言，我相信世界是随机的，很多事情也是不可预测的，我想大部分人也会这样认为，因为我们总是要为自己的存在而辩护，我相信自己在意的那些事物或人有他内在的，无法被计算的珍贵的东西，这种“相信”随着人们对世界和宇宙的不断认识而在不断被加强。
混沌理论的发展也在其中，它是拉普拉斯妖可以被证明为真的路上一座难以逾越的高峰。
故事可以从爱德华·洛伦茨说起（Edward Lorenz），1948年他从麻省理工获得气象学博士学位，之后继续在学校从事气象学研究，1961年末，洛伦茨用学校的计算机去计算自己设计的一组包含12个微分方程的大气模型，当时计算结果保留的精度是小数点后6位。
在同年，他进行了第二次计算，但这次他并没有从头开始，而是选择上次计算时某个中间结果作为初始值，并且只精确到了小数点后3位。
一般人们会认为，同一组方程，差不多精度的初始条件，那么最终得到的结果即使有差别，也应该不会太大，洛伦茨当时可能也是这样认为的吧。然而计算机在经过一段时间之后，他发现输出已经严重偏离了第一次的计算结果，成了两组完全不同的数据。然后他开始研究起这组方程，并把研究结果于1963年在论文上发表出来，文中首次描述了有确定性的动力学系统的混沌现象，从此混沌系统开始进入人们的视野。
简化的洛伦茨方程形式如下：
\[ \begin{aligned} \dot{x} &= \sigma(y - x) \\ \dot{y} &= \rho x - y - xz \\ \dot{z} &= xy - \beta z \end{aligned} \tag{10.1}\]
我们可以把 \(x\)、\(y\) 和 \(z\) 想象成三维坐标下的一个点，\(\dot{x}\)、\(\dot{y}\)、\(\dot{z}\) 分别是它们关于时间的导数， 式 10.1 中的 \(\sigma\)、\(\rho\) 和 \(\beta\) 都为常数且大于 0。
取值 \(\sigma = 10\)，\(\rho = 28\)，\(\beta = \frac{8}{3}\)，针对两组稍微不同的初始条件（[\(x = 1, y = 1, z = 1\)]，[\(x = 1.2, y = 0.8, z = 1.4\)]），分别计算方程组的解，并绘制在三维坐标下。
可以看到两条曲线的轨迹还是有很明显的差异，并且没有收敛（或者固定）于某个特定的值，而是围绕两个中心旋转（也叫不稳定鞍点）：
轨迹在鞍点之间来回游荡，但更重要的是，即使是系数已经确定的微分方程，对于初始条件也可能极其敏感，这种现象被称为蝴蝶效应，蝴蝶扇动翅膀，一个小小的改变可能会带来一场飓风（有趣的是当 \(\rho < 24.74\) 时，两个不稳定点一个会消失，另一个会变成稳定点，轨迹会逐渐收敛到那个稳定点去）。
如果宇宙同样由一组方程描述，那么它应该也会是一个混沌系统，我们不可能完美地无限精度地知道所有初始条件，从这个意义上来说，未来是不可预测的。
但人们也会思考，是不是因为计算机精度有限，才造成这样的结果，于是人们也开始迫切地想要在模拟电路上发现混沌现象（模拟电路有无穷高的精度），虽然后来证明了这不是计算机精度导致的“幻觉”，但人们对混沌电路的研究却也愈发热烈。1983年，任教于加州大学伯克利分校的华裔科学家蔡少棠在日本早稻田大学提出并实现了第一个混沌电路，这个电路被称为“蔡氏电路”，它只由简单的电感、电容、电阻和运放搭建而成，却产生了很明显的混沌效应，暂且不去讨论蔡氏电路，而是转而从基本元素开始，尝试去构建混沌电路。
10.2 混沌电路构建法则
首先这套方法来自《A Concise Guide to Chaotic Electronic Circuits》这本书。
混沌电路其实没有它的名字听起来那么麻烦，我第一次听到混沌电路的时候，觉得很惊讶，它到底是个什么东西。实际上混沌电路只是在模拟微分方程组，有点不好理解？或者说混沌电路本身只是一个计算微分方程结果的计算器。
为了在模拟电路上搭建这个计算器，我们需要先掌握一些基本元素。首先是反相放大器，我们已经很熟悉了，反相放大器反相端可以有多个输入构成加法器：
\[ V_{\text{out}} = - \frac{R_f}{R_1} V_{\text{in}1} - \frac{R_f}{R_2} V_{\text{in}2} - \frac{R_f}{R_3} V_{\text{in}3} \tag{10.2}\]
上式很容易被证明，因为负反馈下反相端处的电压也是 0，如果把同相端输入考虑进去，想要得到类似 式 10.2 的格式，可以考虑这样的结构：
当 \(R_p\) 满足 \(\frac{1}{R_f} + \frac{1}{R_1} + \frac{1}{R_2} = \frac{1}{R_3} + \frac{1}{R_4} + \frac{1}{R_p}\) 时，输入和输出之间满足：
\[ V_{\text{out}} = - \frac{R_f}{R_1} V_{\text{in}1} - \frac{R_f}{R_2} V_{\text{in}2} + \frac{R_f}{R_3} V_{\text{in}3} + \frac{R_f}{R_4} V_{\text{in}4} \tag{10.3}\]
有了加法项之后，还得想办法产生微分项，考虑电容电压和电流是微分关系，对于下面这个电路：
假设运放输出端电压为 \(V\)，有 \(\frac{V - V_{\text{out}}}{R_0} = C \frac{\mathrm{d}V_{\text{out}}}{\mathrm{d}t}\)，两边化简一下可以得到：
\[ C R_0 \dot{V}_{\text{out}} = - V_{\text{out}} + V \tag{10.4}\]
\(V\) 是运放输出端的电压，如果把运放配置为 图 10.2 或 图 10.3 的模式，就可以很好地模拟多输入微分方程了。
如果把电容接在反馈路上，可以构成米勒积分器：
对于高频变化的信号，电容电流满足：
\[ \begin{aligned} \frac{V_{\text{in}}}{R_1} &= - C \frac{\mathrm{d}V_{\text{out}}}{\mathrm{d}t} \\ C R_1 \dot{V}_{\text{out}} &= - V_{\text{in}} \end{aligned} \tag{10.5}\]
有了这些之后，再来处理乘积项，用到的芯片是 AD633，AD633 是一枚乘法器芯片，比起运放来说价格已经非常昂贵了，均价在 20 元左右，除开两个供电引脚外，它有 5 个输入引脚 \(X_1\)、\(X_2\)、\(Y_1\)、\(Y_2\) 和 \(Z\)，以及 1 个输出引脚 \(W\)，满足：
\[ W = \frac{(X_1 - X_2)(Y_1 - Y_2)}{10} + Z \tag{10.6}\]
芯片内部示意图如下：
差不多有了这些东西我们就可以开始搭建电路了，以洛伦茨方程为例，从 图 10.1 的轨迹可以看出，它的 \(x\)、\(y\)、\(z\) 的运动范围比较大，我们第一步是需要把图形整体缩小，不要超出电源电压的范围。
原始方程为：
\[ \begin{aligned} \dot{x} &= 10(y - x) \\ \dot{y} &= 28x - y - xz \\ \dot{z} &= xy - \frac{8}{3}z \end{aligned} \tag{10.7}\]
现在令 \(X = \frac{1}{10}x\)、\(Y = \frac{1}{10}y\)、\(Z = \frac{1}{20}z\)，则有 \(x = 10X\)、\(y = 10Y\)、\(z = 20Z\)，把它们代入到 式 10.7，可以得到：
\[ \begin{aligned} \dot{X} &= 10(Y - X) \\ \dot{Y} &= 28X - Y - 2XZ \\ \dot{Z} &= 5XY - \frac{8}{3}Z \end{aligned} \tag{10.8}\]
这样处理可以使 \(X\)、\(Y\)、\(Z\) 的轨迹与原始 \(x\)、\(y\)、\(z\) 的轨迹结构保持相同，但在每个轴上按照不同比例缩放（也就是线性变换）。
接下来就是用前面介绍的电路分别实现 式 10.8 的微分方程。
对于 \(\dot{X} = 10(Y - X)\)，根据 图 10.4 的电路，改写为 \(\dot{X} = -X + X + 10Y - 10X\)，其中 \(X + 10Y - 10X\) 项等于 式 10.4 的 \(V\)，据此画出电路结构：
此时在电容输出端（有 \(X\) 标记处）的电压满足关系：
\[ 1\text{k}\Omega \cdot 100\text{nF} \cdot \dot{X} = - \frac{10\text{k}}{1\text{k}}X + \frac{10\text{k}}{10\text{k}}X + \frac{10\text{k}}{1\text{k}}Y \tag{10.9}\]
\(1\text{k}\Omega \times 100\text{nF}\) 是一个时间缩放因子，因为 \(\dot{X}\) 的定义是关于时间的微分 \(\dot{X} = \frac{\mathrm{d}X}{\mathrm{d}t}\)，现在定义一个新的“时间” \(\tau\)，令 \(\frac{\mathrm{d}X}{\mathrm{d}\tau} = 1\text{k}\times 100\text{n} \cdot \dot{X}\)，则 \(\tau = \frac{1}{1\text{k}\times 100\text{n}}t\)。
该怎么理解 \(\tau\) 呢，可以把它看成是某种“分辨率”，原本需要花费 \(t\) 时间观测到的动作现在需要花费 \(\frac{1}{1\text{k}\times 100\text{n}}t = \frac{1}{10^{-4}}t = 10^4 t\) 才能观测到（时间延长了），这样我们就能注意到变化过程中更多的细节（由于电子元件自身的特性，过于快速的变化可能会让混沌电路进入稳态，从而失去混沌效应）。
按照这种思路，分别可以把 \(Y\) 和 \(Z\) 的电路搭建出来，然后连接好，上电后自己就能启动了，关于如何搭建以及混沌电路更多细节，可以自行参考开头提到过的那本书。
10.3 更多一点
这一小节不是关于电路的，甚至本章与电路都没太多关系。对于混沌电路，我所了解的也仅限这些了。
但是混沌本身有很多可以说的东西，它被包含在一个大类里面：复杂系统。
对于什么是复杂系统，目前还没有一个很好的定义，不同的学科有它们各自的观点和理解，但就像梅拉妮·米歇尔在她写的书《复杂》里说的，复杂系统都有一个特点，就是在你看到它的第一眼，就知道它是复杂的。
从逻辑斯蒂方程开始吧（Logistic equation），逻辑斯蒂方程最开始是用来描述种群增长规律的（起初是一个微分方程，后来演变成差分方程），假设在某一时刻种群数量为 \(N_t\)，在下一个观测时刻，种群数量 \(N_{t+1}\) 满足：
\[ N_{t+1} = r N_t \left( 1 - \frac{N_t}{K} \right) \tag{10.10}\]
其中 \(r\) 是某比例系数，\(K\) 是环境承载容量，也就是极限值， 式 10.10 表明如果某时刻种群数量达到环境所能承受的最大值，那下一时刻这个种群数量将归零。
对 式 10.10 左右两边分别除以 \(K\)，可以得到一个一般性公式：
\[ x_{n+1} = r x_n (1 - x_n) \tag{10.11}\]
单从这个等式入手，它看起来平平无奇，好像和众多数学公式没有什么差别，也看不出什么奇怪的特性。
现在假设 \(x_0\) 初始值是 0.2，考虑 \(r\) 取不同值的情况，当 \(r = 0\) 时，从 \(x_1\) 开始往后都是 0，随着 \(r\) 逐渐变大，\(x\) 的变化情况如图所示：
在 \(r < 3.8\) 时，\(x_n\) 在迭代一定次数后保持稳定，当 \(r = 3.8\) 时，如果仔细看图的话，会发现图形已经没有规律了，在 \(r = 4\) 时无规律性更明显。
也会有人觉得是因为迭代次数不够，序列还没有收敛？也许它会在第 1000 次或 10000 次以后保持稳定形状？
嗯……但不是这样的，\(r = 4\) 时即使迭代 100 万次，都不会出现稳定或周期性的图案，3.57 是一个分界点，\(r\) 超过 3.57 后整个序列会进入到混沌状态。
就像洛伦茨方程一样，一个完全确定的、看上去平平无奇的方程（式 10.11）也可以出现混沌效应。
现在我们对 \(r\) 连续取值，并把所有 \(r\) 下不断迭代的 \(x\) 具体的值填入以 \(r\) 为横坐标的一条竖线上（相当于统计不同数值出现的次数）：
某个具体数值出现的次数越多，\(Y\) 轴对应那块区域就会颜色越深。把这种映射方式应用到不同的 \(r\)，我们就会得到逻辑斯蒂映射分岔图：
逻辑斯蒂映射分岔图是自相似分形的，意思是我们可以无限放大它，然后会发现放大后的图形与原始图形（图 10.10）具有相似的结构（大概是我是由我自身无穷的叠加形成的）。从图中可以看到，\(r\) 大约在过了 3.6 或 3.7 之后，黑色开始变密，也开始出现混沌效应，如果我们细看 \(r = 3.7\) 处两个更小的分支，会发现它同样也在分岔（随着 \(r\) 的变大， 图 10.8 \(Y\) 轴某些区域会变得更加“混沌”）。
类似的还有帐篷映射（Tent map），它的差分方程为：
\[ x_{n+1} = \mu \cdot \min(x_n, 1 - x_n) \tag{10.12}\]
其中 \(\mu\) 是比例系数，\(\min(a, b)\) 是取最小值运算，这方程看起来比 式 10.11 简单得多，当初始值 \(x_0\) 在 0～1 之间，\(\mu\) 在 0～2 的范围内时，\(x\) 总可以保持在 0～1 的范围内，然而当 \(\mu > 1.4\) 时，输出序列会进入到混沌状态：
