自然的語(yǔ)言究竟是什么？

思考一下這個(gè)問(wèn)題。如果你問(wèn)一個(gè)物理學(xué)家，我敢打賭他們的回答會(huì)是微分方程（differential equation）的某種變體。自然界中的一切，從原子的振動(dòng)，到颶風(fēng)的形成，再到行星和恒星在浩瀚太空中的運(yùn)動(dòng)，都是由微分方程控制的。

實(shí)際上，流體的混沌流動(dòng)是由納維-斯托克斯方程（Navier-Stokes）微分方程支配的，甚至牛頓的第二定律F=ma也是一個(gè)微分方程，因?yàn)榧铀俣仁俏恢玫亩A導(dǎo)數(shù)。

這個(gè)想法宏大而美麗，我們可以用一個(gè)相對(duì)簡(jiǎn)單的數(shù)學(xué)概念來(lái)描述周圍的許多事物。然而，有一個(gè)不小的問(wèn)題：大多數(shù)微分方程根本無(wú)法求解。自然界的所有秘密都被它們的數(shù)學(xué)復(fù)雜性完全封鎖住了，或者說(shuō)真的如此嗎？雖然大多數(shù)方程不可解，但可以通過(guò)迭代找到近似解。這個(gè)主題非常廣泛且復(fù)雜。不過(guò)現(xiàn)在我們先聚焦一個(gè)簡(jiǎn)單的例子。

一個(gè)在彈簧上來(lái)回運(yùn)動(dòng)的方塊，這是一個(gè)經(jīng)典的簡(jiǎn)諧振子（Harmonic Oscillator）例子。事實(shí)上，它足夠簡(jiǎn)單，以至于描述該系統(tǒng)的微分方程確實(shí)有一個(gè)確切的解。這其實(shí)就是牛頓的第二定律——力等于質(zhì)量乘以加速度，再結(jié)合彈簧力定律——力等于負(fù)的彈簧常數(shù)乘以位移。

通過(guò)重新排列，可以解出加速度。當(dāng)解這個(gè)微分方程時(shí)，會(huì)得到一個(gè)關(guān)于方塊位移的正弦波解，隨著時(shí)間變化而變化。

由于位移是正弦波形式的，所以速度也是。如果我們知道所有的參數(shù)，比如振幅、相位因子（即啟動(dòng)時(shí)鐘的位置）、方塊的質(zhì)量和彈簧的剛度，那么就可以精確預(yù)測(cè)在任何時(shí)刻方塊的位置。

不過(guò)假設(shè)你不知道如何解牛頓方程，那你該如何近似模擬這個(gè)系統(tǒng)的運(yùn)動(dòng)呢？通過(guò)使用牛頓第二定律和胡克定律，我們發(fā)現(xiàn)系統(tǒng)的加速度為負(fù)的彈簧常數(shù)除以質(zhì)量再乘以位移。

如果將連續(xù)時(shí)間分割成離散的時(shí)間步長(zhǎng)ΔT，就可以通過(guò)迭代計(jì)算速度v和位移x。

在每一個(gè)時(shí)間步長(zhǎng)上，力會(huì)增加或減少方塊的速度，

速度也會(huì)相應(yīng)地增量位置,

每個(gè)時(shí)間步長(zhǎng)，位移、速度和力都會(huì)發(fā)生變化，因此每次都需要重新計(jì)算。

讓我們用這個(gè)方法來(lái)近似位移與時(shí)間的正弦曲線，

紅線是精確解。你可以看到一開始近似還是挺好的，但隨著時(shí)間推移，準(zhǔn)確性迅速降低。這是任何近似方法的常見問(wèn)題。然而還有其他的問(wèn)題。注意到近似的振幅隨著時(shí)間增加在增大，盡管系統(tǒng)是封閉的，這意味著彈簧在每個(gè)周期都會(huì)越來(lái)越伸展，顯然這是不合邏輯的，而且也違反了能量守恒。

這其實(shí)是我們所用方法的一個(gè)特定缺陷——它不太能很好地守恒量，比如能量。經(jīng)過(guò)多次迭代后，近似結(jié)果可能會(huì)趨向無(wú)限大。

通過(guò)選擇較小的時(shí)間步長(zhǎng)，可以緩解這個(gè)問(wèn)題。步長(zhǎng)越小，誤差積累得越慢。不過(guò)，這也有問(wèn)題，較小的步長(zhǎng)會(huì)增加計(jì)算量，因?yàn)樾枰?jì)算更多的迭代步驟。另一方面，能有一種能處理較大時(shí)間步長(zhǎng)的技術(shù)會(huì)更好，這樣既提高計(jì)算效率，也能在不需要過(guò)多小步長(zhǎng)的情況下推進(jìn)系統(tǒng)的演算。這種近似方法有其局限性，但通過(guò)調(diào)整時(shí)間步長(zhǎng)和方法，我們可以逐步改善結(jié)果，并更好地模擬自然現(xiàn)象。

能否找到一種沒(méi)有當(dāng)前方法問(wèn)題的方式？讓我們用更通用的形式寫出隱式歐拉方法。

首先，我們實(shí)際上是在解一個(gè)由兩個(gè)微分方程組成的系統(tǒng)，但為了簡(jiǎn)化，現(xiàn)在我們只討論一個(gè)一般性的方程，我們?cè)噲D求解的方程形式是

帶有某些初始條件

之前使用的方法實(shí)際上有個(gè)名字，叫歐拉法 (Euler's Method)。由于有初始條件，我們通常稱這種問(wèn)題為初值問(wèn)題 (Initial Value Problems，簡(jiǎn)稱 IVP)。

再一次，如果我們解析地求解這個(gè)方程，那么就將時(shí)間t離散化為由小時(shí)間步長(zhǎng) Δt的集合。然后，通過(guò)歐拉法，如果知道某些初始條件x_i在時(shí)間t_i，那么通過(guò)微分方程，可以找到x_(i+1)的近似值

之所以能這么做，是因?yàn)楦鶕?jù)原方程，函數(shù) f(x,t)是該點(diǎn)處的導(dǎo)數(shù)，

因此我們本質(zhì)上是在做一個(gè)微小的線性近似。

如前所述，時(shí)間步長(zhǎng)足夠小的情況下，這種方法有效。然而如果時(shí)間步長(zhǎng)過(guò)大，會(huì)遇到問(wèn)題。我們可以對(duì)歐拉法做一些小的調(diào)整來(lái)解決這些問(wèn)題。通過(guò)隱式歐拉法，我們?cè)谟?jì)算x_(i+1)時(shí)，利用的是

而不是

通過(guò)這種方式，或許可以減少之前遇到的問(wèn)題，因?yàn)槲覀冊(cè)谟?jì)算中考慮了下一個(gè)位置。

注意，這種數(shù)值方法將更加數(shù)值穩(wěn)定。我們稱這種方法為隱式的，因?yàn)閤_(i+1)依賴于自身。而之前介紹的歐拉法可以被視為顯式方法。

GIF中從上到下顯示的是：精確解、顯式歐拉法和隱式歐拉法。為什么隱式彈簧慢慢停止了？因?yàn)槲覀冊(cè)谟?jì)算中沒(méi)有包含阻尼項(xiàng)，這就是我們做出的取舍。隱式方法比顯式方法數(shù)值上更穩(wěn)定，但換來(lái)的代價(jià)是每個(gè)時(shí)間步長(zhǎng)會(huì)漏出一些能量，導(dǎo)致系統(tǒng)表現(xiàn)得像是有阻尼一樣。實(shí)際上，如果我們重新排列隱式方法，就會(huì)發(fā)現(xiàn)它與顯式方法的計(jì)算是反過(guò)來(lái)的。我們從x_(i+1)計(jì)算x_i，因此在某種意義上，隱式方法就是一個(gè)反向的歐拉法。從這個(gè)角度理解，就很容易明白為什么會(huì)有阻尼的現(xiàn)象。

顯式歐拉法隨著時(shí)間推移振幅會(huì)增大，而隱式方法則相反，振幅隨著時(shí)間減少。因此，這種近似方法并不適合彈簧模擬。

那么我們?nèi)绾谓Y(jié)合顯式和隱式方法的優(yōu)點(diǎn)，同時(shí)弱化它們的缺點(diǎn)呢？最簡(jiǎn)單和最直接的方式可能是取它們的平均值，雖然不完全是，但想法很類似。我們這么做的目的是希望顯式方法的能量增加可以抵消隱式方法的能量減少。為此，我們計(jì)算兩個(gè)值K_1和K_2，

注意我們?cè)谟?jì)算K_2時(shí)使用的是顯式歐拉法中的結(jié)果。然后我們將這兩個(gè)值的平均值加到x_i 上，這就得到了更精確的結(jié)果。

這一方法能在保持?jǐn)?shù)值穩(wěn)定性的同時(shí)，也較好地守恒能量。

這就是我們的方法了，它比之前的方法更加數(shù)值穩(wěn)定，同時(shí)能量守恒也更好。看一看，幾乎完全吻合，

不過(guò)我們可以看出這仍然是一個(gè)近似解，因?yàn)殡S著時(shí)間的推移，準(zhǔn)確度還是會(huì)下降，但這無(wú)疑是目前最好的方法。實(shí)際上，我們剛剛使用的這個(gè)技巧有一個(gè)專門的名稱，叫做二階龍格-庫(kù)塔方法 (Runge-Kutta 2nd Order，RK2)。這個(gè)方法被稱為“二階”是因?yàn)槲覀兤骄藘蓚€(gè)項(xiàng)K_1和K_2。事實(shí)上，顯式歐拉法等同于一階龍格-庫(kù)塔方法，這就是為什么顯式歐拉法中的公式也會(huì)出現(xiàn)在我們對(duì) RK2 的計(jì)算中。

同樣需要注意的是，龍格-庫(kù)塔方法既可以是顯式的，也可以是隱式的?，F(xiàn)在，實(shí)際應(yīng)用中最常用的是四階龍格-庫(kù)塔方法，簡(jiǎn)稱 RK4，它也是基于相同的概念，只是加了幾個(gè)額外的權(quán)重或 K 值。

在這里，我們不僅關(guān)心t_i和 t(i+1)，還關(guān)心它們之間的中點(diǎn)。和之前一樣，我們計(jì)算K_1，并將其加到x_i上，然后將其外推到中間點(diǎn)。在該點(diǎn)計(jì)算K_2，并再次用該斜率找到另一個(gè)中間點(diǎn)，接著計(jì)算K_3。最后，用K_3斜率完成一個(gè)完整的時(shí)間步長(zhǎng)，并在該點(diǎn)計(jì)算K_4。

然后平均這四個(gè)斜率值，再乘以 Δt，得到 x(i+1)。

這里是完整的代數(shù)公式，

注意K_2和K_3的權(quán)重比K_1和K_4更高。這種權(quán)重分配來(lái)源于 RK4 的推導(dǎo)過(guò)程，其中涉及到泰勒級(jí)數(shù)。我更多的是想建立一種直觀理解，說(shuō)明為什么以及如何使用龍格-庫(kù)塔方法。

這就是我們最終構(gòu)建的近似器，RK4 廣泛應(yīng)用于科學(xué)和工程領(lǐng)域的初值問(wèn)題。其他龍格-庫(kù)塔方法也被廣泛使用，例如Matlab的ODE45 就使用了四階和五階龍格-庫(kù)塔方法的組合。你會(huì)注意到，雖然 RK2 方法已經(jīng)足夠好了，但 RK4 方法更加精確，隨著時(shí)間推移，兩者的準(zhǔn)確度差距會(huì)越來(lái)越明顯。具體使用哪種近似方法取決于應(yīng)用場(chǎng)景，對(duì)于某些系統(tǒng)，顯式或隱式歐拉法就足夠了，而 RK4 通常是首選，因?yàn)樗諗克俣瓤烨铱傮w上更穩(wěn)定。

求解初值問(wèn)題的近似解是數(shù)值分析領(lǐng)域中一個(gè)極為廣泛的課題。但僅使用我概述的這些原理，你已經(jīng)可以創(chuàng)造許多很酷的東西了。