最初計算機出現就是為解決兩方面問題：第一破解密碼，第二科學計算。所以想聊聊這兩方面的內容，也就是數學、物理學。有時候換個角度會發現自己的技能樹還能解決其他學科的問題。而這些都是很實用的問題，其衍生出了衛星定軌、天氣預報、礦產勘探等一系列內容。

本意想模仿《時間簡史》的筆法，但力有未逮。所以次分享除了科普，還有實現部分，甚至於還有公式（少量解釋問題使用）。

本次主要聊聊非線性系統。混沌問題實際上是一個非線性方程，它可以用形象的諸如“蝴蝶暴雨”之類的比喻解釋，也可以用“機器學習”研究其統計特徵。因此實現部分首先用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 [email protected]')    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）：

目前在水流、冰雪模擬方面走在前沿的是一些動畫公司，比如獲得奧斯卡最佳動畫短篇的《鷸》：

感興趣的可以體驗一下迪士尼的炫技。可以想象每一幀的背後需要多少計算資源。

對於連續的非線性系統模擬有幾個比較常用的思想：第一是用小體積近似思想，由此產生出了有限元等方法；第二是使用粒子近似，比較有代表性的是 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 會員，暢享 1000+ 場 Chat ！點選檢視