写给程序员的数理科普：混沌与三体

最初计算机出现就是为解决两方面问题：第一破解密码，第二科学计算。所以想聊聊这两方面的内容，也就是数学、物理学。有时候换个角度会发现自己的技能树还能解决其他学科的问题。而这些都是很实用的问题，其衍生出了卫星定轨、天气预报、矿产勘探等一系列内容。

本意想模仿《时间简史》的笔法，但力有未逮。所以次分享除了科普，还有实现部分，甚至于还有公式（少量解释问题使用）。

本次主要聊聊非线性系统。混沌问题实际上是一个非线性方程，它可以用形象的诸如“蝴蝶暴雨”之类的比喻解释，也可以用“机器学习”研究其统计特征。因此实现部分首先用TensorFlow的RNN函数完成了一个LC电路模拟，还有一个三体问题模拟，之后用“混沌吸引子”绘制几张壁纸。

1. 前言：从 LC 震荡电路到 RNN 神经网络

几乎每篇文章都会写前言，在前言里会聊聊对于文章内容的看法。本次文章的内容集中于数理科普，依然有机器学习的内容，因为机器学习与数值模拟总是不分家的或者说所有的自然科学都是有迹可循的。了解一下其他学科内容权当开阔思路。

以前听高晓松炫耀过，他们家里人如何如何的聪明，自己焊个收音机啥的。其实这并不难，而且本应该是理工科的天赋技能，并没有什么值得炫耀的。虽然这其中会涉及模电数电之类的知识，比如电容、比如三极管，但是并不是不可逾越的障碍。所谓了解意味着知道每种电器原件对于电流和电压的响应。使用最简单的 LC 震荡电路来讲，其大概的形状是这样的：

这里就有可说的了，高中就学过，LC 震荡电路实际上可以产生正弦波，对于理想情况（没有内阻）可以用一个二次常微分方程约束：

$$\begin{matrix}i=C\frac{du}{dt}\u=L\frac{di}{dt}\end{matrix}\rightarrow u=-LC \frac{d^2u}{dt^2}$$

上式中 $i$ 是电流、$u$ 是电压、$C$ 是电容大小、$L$ 是电感大小。关于电压 u 的常微分方程的解实际上就是一个三角函数。在数值求解的过程中，一般会这么做：

$$\begin{matrix}\dot{u}=y\\dot{y}=m\cdot u\\frac{df}{dt}\rightarrow\dot{f}\m\rightarrow- \frac{1}{LC}\end{matrix}\rightarrow \begin{matrix}x\rightarrow [u, y]^T\\dot{x}=A\cdot x\end{matrix}$$

上面的公式利用了几个思想：第一个降低方程次数；第二个转化为矩阵表示。看到这种时序依赖问题就会想到 RNN，但是直到现在还没有看出来这个公式和 RNN 有什么联系。下面需要很重要的一步，就是“差分”。有限差分是数值模拟里用的最多的一种方法，其使用了近似的思想：

$$\dot{x}\approx \frac{x_{t+1}-x_{t-1}}{2dt}\approx A\cdot x_t$$

稍微整理一下

$$x_{t+1}=M\cdot [x_t, x_{t-1}]^T$$

这几乎就是一个不带激活函数的循环神经网络。因此可以使用 RNN 网络去完成 LC 电路的数值模拟工作。将上述公式整理成与 RNN 相同的形式就可以进行模拟了，这是使用 TensorFlow 完成的 LC 震荡电路模拟：

import numpy as npimport tensorflow as tffrom tensorflow.contrib.rnn import RNNCellclass MyRNN(RNNCell):    """    定义无激活函数的RNN    """    def __init__(self, w):        super(MyRNN, self).__init__()        self.W = w    def call(self, inputs, state):        x = tf.concat([inputs, state], 1)        ret = tf.matmul(x, self.W)        return ret, retLC = 10dt = 0.01M = np.array([[0, 1 * dt], [-1/LC * dt, 0], [1, 0], [0, 1]])graph = tf.Graph()with graph.as_default():    inputs = tf.placeholder(tf.float32, [1, 2])    states = tf.placeholder(tf.float32, [1, 2])    W = tf.Variable(M, dtype=tf.float32)    cell = MyRNN(W)    output, state = cell(inputs, states)    init = tf.global_variables_initializer()sess = tf.Session(graph=graph)sess.run(init)xt = np.array([[0, 2]])xtn1 = np.array([[0, 2]])lines = np.zeros([20000])t = []for itr in range(20000):    xtp1, xtp1 = sess.run([output, state], feed_dict={inputs:xt, states:xtn1})    xtn1 = xt    xt = xtp1    lines[itr] = xt[0, 0]    t.append(dt/2 * itr)import matplotlib.pyplot as pltimport matplotlibmatplotlib.style.use('ggplot')plt.rcParams['font.sans-serif'] = ['Source Han Serif CN']plt.rcParams['axes.unicode_minus'] = Falseplt.text(0,-1.5,u"@如是")plt.plot(t, lines, lw=2, alpha=0.5, color="#990000")plt.show()

来看下模拟效果：

其实 LC 电路电压、电流均是三角函数，其周期和 LC 大小有关。震荡电路经常用于钟表之中，收音机天线中也用到了。我们的 RNN 还有一个激活函数，其形式如下，它带来了非线性：

$$\dot{x}=f(A\cdot x)$$

2 引入非线性-混沌吸引子

2.1 混沌问题引入

首先来看最出名的洛伦兹吸引子图形：

其方程形式：

$$\begin{matrix}\dot{x}=10(y-x)\\dot{y}=x(28-z)\\dot{z}=xy-\frac{8}{3}z\\end{matrix}\rightarrow \dot{\vec{x}}=f(\vec{x})$$

可以看到，洛伦兹吸引子的方程相当于在前面的线性微分方程的基础上添加了非线性项。因此得到了一个看起来很奇怪的图形。这种结果有的人称之为“混沌”，混沌系统对于数值精度和初始值都是敏感的。稍微改变一下位置初始值：

就会得到不太不同的结果。代码：

"""洛伦兹吸引子"""from mayavi import mlabimport numpy as npfrom scipy.integrate import odeint# 这里odeint就是上面RNN代码细化def lorenz(w, t, p, r, b):    x, y, z = w    return np.array([p*(y-x), x*(r-z)-y, x*y-b*z])t = np.arange(0,36, 0.01)track1 = odeint(lorenz, (0., 15., 0.), t, args=(10., 28., 3.))track2 = odeint(lorenz, (0., 10, 0.1), t, args=(10., 28., 3.))mlab.plot3d(track1[:,0], track1[:,1], track1[:,2], color=(0.7,0,0), tube_radius=0.5)mlab.plot3d(track2[:,0], track2[:,1], track2[:,2], color=(0.0,0.7,0), tube_radius=0.5)mlab.show()

如何解释这个问题呢：如果我们站在三维空间的某一点，这个点上有一个小的箭头，你沿着箭头走一步，接下来再沿着当前箭头走一步，那么最终你走过的轨迹就是我们所谓的常微分方程的解。而空间中一系列的为我们所指定方向的箭头称之为“场”。对于洛伦兹吸引子其场的形式如下：

可见场的形式非常复杂，这意味着我们走偏差一些就会使得最终的终点完全不同。对于线性常微分方程则不具有这种敏感性：

线性问题的方程其实就是 LC 电路那个形式：

$$\dot{\vec{x}}=A\cdot \vec{x}$$

当然，还有其他形式的吸引子，比如：

其数学形式为：

$$\begin{matrix}\dot{x}=40(y-x)+0.16xz\\dot{y}=55x-xz+20y\\dot{z}=1.833z+xy-0.65x^2\end{matrix}$$

代码：

"""吸引子"""from mayavi import mlabimport numpy as npfrom scipy.integrate import odeintdef attrack(w,t,a1,a2,a3,a4,a5,a6):     x, y, z = w    x = x    return np.array([a1*(y-x)+a4*x*z,a2*x-x*z+a6*y,a3*z+x*y-a5*x*x])t = np.arange(0, 30, 0.0001)track = []track1 = (odeint(attrack, (1.1, 0, 0.0), t, args=(40,55,1.833,0.16,0.65,20))/300+0.5)mlab.plot3d(track1[:,0], track1[:,1], track1[:,2], color=(0.7,0.,0.), tube_radius=0.01)mlab.show()

这个代码调整一下曲线的透明度和粒子显示就可以得到上面图形的效果。

2.2 蝴蝶煽断翅膀也引不来飓风

混沌系统对于初始值和计算精度十分敏感，这使得一个小的扰动会很大的改变系统的最终状态。于是那个蝴蝶就来了：

亚马逊雨林一只蝴蝶翅膀偶尔振动，也许两周后就会引起美国得克萨斯州的一场龙卷风。

本身而言混沌系统难于精确模拟，但其是有规律可循的。系统本身最终也会趋于一个稳态，这种稳态也就是吸引子。从另一个角度来说，这是混沌系统本身的统计特征。统计的引入使得我们可以使用机器学习算法去研究。

混沌系统虽然难以得到精确解，但其并不太可能会出现不合常理的结果。因此虽然最终状态随着初始值会有所不同，但最终的状态是可以通过统计来预知的。“蝴蝶效应”本身可能并没有描述的那么恐怖，最大的可能就是蝴蝶煽动了一下翅膀，然后掉了一片树叶，而这片树叶是被翅膀拍下去的。

3. 卫星与三体问题

3.1 卫星轨道

中学时学过的万有引力定律实际上就是一个非常简单的常微分方程。在求解卫星轨道的过程中进行了相当程度的简化，比如地球是不动的，卫星仅受到地球引力的作用，考虑到人造卫星非常小而且离地球非常近，这种简化通常是合理的（如果要没有月球的话）。这种约束与我们上面提到的 LC 震荡电路十分类似，或者说这根本上就是一个二次常微分方程：

$$\ddot{\vec{x}}=\frac{GM}{\vec{x}^2}\hat{e}_r$$

使用同样的思路降次求解。可以得到：

求解代码：

"""轨道问题求解"""import numpy as npfrom scipy.integrate import odeintimport matplotlib.pyplot as pltimport matplotlibmatplotlib.style.use("ggplot")def orib(w, t, GM):    # x1, x2 坐标    # v1, v2 速度    x1, x2, v1, v2 = w    sqrt3 = np.sqrt(x1 ** 2 + x2 ** 2)    # dx / dt = v    x1n = v1    x2n = v2    # dv / dt = GM/r^2    v1n = - GM * x1 / sqrt3    v2n = - GM * x2 / sqrt3    return [x1n, x2n, v1n, v2n]# 给定时间t = np.arange(0, 300, 0.001)# (1, 2, 3, 4)初始位置和速度data = odeint(orib, (1, 2, 3, 4), t, args=(5, ))# 画图plt.plot(data[:, 0], data[:, 1], lw=2, alpha=0.5)plt.axis("equal")plt.show()

3.2 三体问题

假设有三个物体仅在万有引力的作用下运动，那么某一时刻以单独的星体来看，其会受到其他两个物体的作用，他们的运动轨迹无法有一个解析的表达，但是可以用常微分方程表示。跪拜一下牛顿，历史上最伟大的科学家没有之一，有兴趣的一定要看看《自然哲学的数学原理》。下面进行一个简单的三体问题模拟：

图片比较大可能加载比较慢，我们稍微改变一下右上角的点的位置，从 2 加到 2.1 来看下结果：

运动情况就完全不同了，这是最简单也是最常提起的“三体问题”，三体作为一个非常常见的“混沌问题”，其对于初始值选取是十分敏感的，也就是我们选择初始值稍微改变，就会使得其后的运动产生很大变化。从图形里可以看到，三个星体随时可能撞到一起，这是很难通过计算预知的。

求解与绘图程序：

# -*- coding: utf-8 -*-"""三体问题求解"""from scipy.integrate import odeint import numpy as np from moviepy.video.io.bindings import mplfig_to_npimageimport moviepy.editor as mpydef Three(w,t,M1,M2,M3):     """    三体问题微分方程    """    # 三个星体的速度和位置    v1x,v1y,v2x,v2y,v3x,v3y,x1,y1,x2,y2,x3,y3= w    # 星体距离    r12=np.sqrt((x1-x2)**2+(y1-y2)**2)    r23=np.sqrt((x3-x2)**2+(y3-y2)**2)    r13=np.sqrt((x1-x3)**2+(y1-y3)**2)    # 万有引力定律    rv1x=M2*(x2-x1)/r12**3+M3*(x3-x1)/r13**3    rv1y=M2*(y2-y1)/r12**3+M3*(y3-y1)/r13**3    rv2x=M1*(x1-x2)/r12**3+M3*(x3-x2)/r23**3    rv2y=M1*(y1-y2)/r12**3+M3*(y3-y2)/r23**3    rv3x=M2*(x2-x3)/r23**3+M1*(x1-x3)/r13**3    rv3y=M2*(y2-y3)/r23**3+M1*(y1-y3)/r13**3    # 速度与位置关系    rx1=v1x    ry1=v1y    rx2=v2x    ry2=v2y    rx3=v3x    ry3=v3y    return np.array([rv1x,rv1y,rv2x,rv2y,rv3x,rv3y,rx1,ry1,rx2,ry2,rx3,ry3]) t = np.arange(0,6000, 0.01) # 创建时间点 # 调用ode对三体问题进行求解track1 = odeint(Three, (0,-2,0.1,0.2,0.21,0.1,2,2.1,-2,-2,-2,2), t, args=(6.4,6.4,7.4)) # 绘图import matplotlib.pyplot as plt import matplotlibmatplotlib.style.use("ggplot")fig = plt.figure()ax=fig.add_subplot(111);ax.set_xlabel("X(km)")ax.set_ylabel("Y(km)")ax.set_title("Three-body")plt.axis('equal')def update_lines(dt):    plt.cla()    d=int(dt*600)    #plt.xlim([-50,50])    #plt.ylim([-50,50])    plt.text(4,6, 'T='+str(dt)+'s')    plt.text(-2,-6, r'Three-body Problem@RS')    ax.plot(track1[d:d+1, 6], track1[d:d+1, 7],'o')    ax.plot(track1[d:d+1, 8], track1[d:d+1, 9],'o')    ax.plot(track1[d:d+1, 10], track1[d:d+1, 11],'o')        ax.plot(track1[0:d, 6], track1[0:d, 7],color='cornflowerblue',lw=2,alpha=0.5)    ax.plot(track1[0:d, 8], track1[0:d, 9],color='orange',lw=2,alpha=0.5)    ax.plot(track1[0:d, 10], track1[0:d, 11],color='b',lw=2,alpha=0.5)    return mplfig_to_npimage(fig)# 输出动画animation =mpy.VideoClip(update_lines, duration=39)animation.write_gif("long.gif", fps=5)

这种对于初始值的敏感性使得我们求解“一个”数值解的过程中会产生问题，因为计算机精度并不是无限精度的。如果是无限精度的话意味着，我们将所有的初始条件考虑进去，是可以预知最终状态的。但是混沌问题的敏感性加之计算机有限的计算精度使得我们几乎无法预知最终状态。这就陷入了“可知论”与“不可知”论则争吵。理论上我们知道事物的开始是可以知道最终结果的，但是对于复杂系统，由于初始条件的复杂性，使得我们几乎无法可以完全预知未来。

3.3 更复杂的自组织临界系统

“蝴蝶效应”所描述的情况，更加类似于“自组织临界系统”。所谓自组织临界系的概念在 1987 年由贝克、唐和威圣菲尔德首次提出来，它用来描述状态过渡，比如，在除了过渡温度之外的所有温度上，系统的扰动仅局部地影响系统的成分。在临界温度，即使只是最邻近的系统成分直接地相互作用，扰动也会影响整个系统。自组织临界系统在群体智能、社会学等都有应用。

举个例子来说，一些研究人员认为菲律宾海啸与汶川地震有所联系。这是因为地下环境的复杂性使得地层、地下水、地热综合起来形成了一个自组织临界系统，而这个系统是全球性的，某个地区发生火山、地震可能会引起地球任何一个地方的地震。如果真是这样的话，那么地震预报基本上是不可能完成的任务。因为你永远无法知道，哪根是压倒骆驼的最后一根稻草。

4. 非线性系统其他应用

在进行数值模拟的时候，非常不喜欢做的就是对于非线性系统的模拟，因为它非常难。

与我们生活关系最大的非线性系统就是流体。有一种说法是目前这个时间节点，传统物理领域唯一可能产生诺贝尔奖的研究就是流体力学领域的。有多么难呢？请看下面的图片：

上面这章图片节选自特仑苏的广告，这叫做“牛奶皇冠”。实际上，指导目前为止都没有一个非常好的算法来形成上面的图形。同时流体模拟需要非常高的计算资源，因此很久以前，你很难在动画中看到水流，或者看到了也是类似于油的质感看下图（参考文献 1）：

目前在水流、冰雪模拟方面走在前沿的是一些动画公司，比如获得奥斯卡最佳动画短篇的《鹬》：

鹬 Piper (2016)【奥斯卡最佳动画短篇】

感兴趣的可以体验一下迪士尼的炫技。可以想象每一帧的背后需要多少计算资源。

对于连续的非线性系统模拟有几个比较常用的思想：第一是用小体积近似思想，由此产生出了有限元等方法；第二是使用粒子近似，比较有代表性的是 SPH、MPM 方法。现在在一些冰雪中使用 MPM 方法比较多。来看看《疯狂动物城》的场景：

雪地下陷的场景就是使用的 MPM 方法进行的模拟。当然这种方法在《冰雪奇缘》中用的更多。可能有些人对计算量级没有概念，我们将物体看成是一个个的小球，小球带有几个参数（浮点数），那么进行最基本的计算需要百万量级以上的粒子。然后再迭代计算各个粒子之间的交互，想想都是一个十分耗时的过程。看下面的图形：

为了模拟这个雪崩，我使用了大概 100×100*50 个粒子（SPH方法），计算了很久。而动画中需要的量级远比这个大的多。而这只是一个简单的连续介质的模拟。所以说动画里看的不是水，而是大把大把的钱。

参考文献

1. Patkar S, Aanjaneya M, Karpman D, et al. A hybrid Lagrangian-Eulerian formulation for bubble generation and dynamics[C]// ACM Siggraph/eurographics Symposium on Computer Animation. ACM, 2013:105-114.

---

本文首发于GitChat，未经授权不得转载，转载需与GitChat联系。

阅读全文: http://gitbook.cn/gitchat/activity/5b28c7f63130202ab558cb48

您还可以下载 CSDN 旗下精品原创内容社区 GitChat App ，阅读更多 GitChat 专享技术内容哦。