在學(xué)習(xí)高數(shù)的過程中，相信很多人都“清楚”地認(rèn)識到大部分積分式不可求解 （準(zhǔn)確來說是懶得解），積分式內(nèi)的結(jié)構(gòu)就是一個dy/dx = f(x)的最基本的微分方程。

但是現(xiàn)實(shí)是大部分一階微分方程都是dy/dx?= f(x, y)的結(jié)果，也就是說目標(biāo)方程本身就在其導(dǎo)數(shù)內(nèi)，such：dy/dx=y，自覺上(查表)知道y=αe^t，但是如果真的要求解的話也會讓很多人無從下手。雖然有各種求解微分方程的辦法，但是都異常繁瑣且不通用。這時候就到了數(shù)值積分發(fā)揮用處的時候了。

什么是數(shù)值積分

為了更貼合使用數(shù)值積分的情況，我們來想象一個情景：一輛汽車在直路上行走，路上有一個原點(diǎn)，并且汽車的速度與汽車到原點(diǎn)的距離成正比，在開始時(t=0)，汽車距離原點(diǎn)1米，求汽車與原點(diǎn)的距離y與時間t的關(guān)系。

根據(jù)題意不難列出：v=dy/dt=αy, y|t=0 = 1，當(dāng)然也不難知道y(t)=e^αt，但是假如我們特別蠢，完全不會求解這個方程，但是又特別想知道結(jié)果是怎么樣的，那么就只能暴力計(jì)算：算出當(dāng)前速度，把汽車按照這個速度前進(jìn)一點(diǎn)點(diǎn)距離，然后再計(jì)算速度，然后再前進(jìn)一點(diǎn)點(diǎn)，然后......

這種做法就是數(shù)值積分的核心思想了，數(shù)值積分的核心思路就是積分的原理:

這時候就會有人說，每一步都會累積誤差，完全不是精確的，怎么能叫求解呢。問得好，盡管是有誤差，但是精度可以控制在接受范圍內(nèi)，科學(xué)家們已經(jīng)用這種辦法對幾百億個星體進(jìn)行演算模擬出星系碰撞甚至宇宙的演化。

了解了數(shù)值積分是什么后，就可以開始了解求解細(xì)節(jié)和優(yōu)化的辦法了

基礎(chǔ)? --? ?歐拉法

雖然第一個是歐拉發(fā)明的方法，但卻是最簡單的一個辦法

注意：以下所有方程都為? dy/dt=f(t,y)??的結(jié)構(gòu)，且附帶初值條件y0=○。目標(biāo)方程y(t)可以為向量函數(shù)，也就是說所有方法都適用于高維空間

看到導(dǎo)數(shù)的定義式：

可以得到：

我們記y_n為當(dāng)前知道的點(diǎn)，對應(yīng)的時間為t_n，y_n+1為下一個要求的點(diǎn)，那么算法可以寫為：

進(jìn)階? --? 中點(diǎn)法

有些人可能看出來了一個問題：如果遇到曲率很大的地方，直接采用y_n點(diǎn)的導(dǎo)數(shù)會與實(shí)際曲線存在很大偏差

對此，我們又想到了一樣?xùn)|西：一段圓弧的中點(diǎn)的切向量總與兩端點(diǎn)連線的向量平行

于是得到中點(diǎn)法方法:

對于歐拉法和中點(diǎn)法分別對上面汽車的例子作圖來看看

可以中點(diǎn)法比歐拉法更接近目標(biāo)方程, 但是也有比較大的誤差.

接下來的方法已經(jīng)沒有明顯的對比效果了，在最后我會用“雙蝴蝶”做一個總結(jié)

究極? --? Heun法（梯形法）

看到不知道什么什么定理（我是真的忘了名字），F(xiàn)是f積分后的原函數(shù)

因?yàn)椴荒艽_定ξ的值，所以一般用ab兩點(diǎn)f的平均數(shù)代替f(ξ)，即：

根據(jù)這個思路得到：

這個方法對一種叫做剛性系統(tǒng)的東西有特攻，后面再詳細(xì)講述, 并且善于利用的話可以做到任意小的誤差

我也不知道是什么原理的方法? --? Runge-Kutta法（rk4）

rk4的計(jì)算方法如下：k1是y_n起點(diǎn)的導(dǎo)數(shù)，k2的第一中點(diǎn)的導(dǎo)數(shù)，k3是第二中點(diǎn)的導(dǎo)數(shù)，k4是k3終點(diǎn)的導(dǎo)數(shù)，各自加權(quán)后得到最終斜率

歐拉法的誤差為Δt^2，梯形法的誤差為Δt^3，而rk4的誤差達(dá)到了Δt^5

因?yàn)橛?jì)算量少，并且誤差小，所以rk4是一棟非常常用的求解微分方程的方法

Butcher tableau? ?和? ?各種顯式rk方法

受rk4的啟發(fā)，人們開始拓展多點(diǎn)權(quán)重累加的方法，通式如下：

為了簡潔描述，發(fā)明了一種叫Butcher tableau的表格：

有了這種表格可以自己嘗試調(diào)整系數(shù), 去尋找快捷方便的求解方法. 但是為了保證方法是自洽的(也就是不會自身出問題), 需要有以下條件:

簡單來說就是：每一行的a相加等于c，且b相加等于1

以下為各種不同的顯式rk方法? en.m.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods

上述方法在大部分情況都表現(xiàn)不錯，但是缺點(diǎn)是需要手動調(diào)整時間間隔以確保性能精度兼得。以懶著稱的數(shù)學(xué)家發(fā)明了根據(jù)誤差自動調(diào)整間隔的自適應(yīng)方法

誤差e的計(jì)算就是用不同的權(quán)重計(jì)算y_n+1，如果兩個結(jié)果相差太大則減少間隔以確保精度。

這里我需要說一聲抱歉，我一直都找不到如何根據(jù)誤差調(diào)整間隔，如果有讀者知道的話請務(wù)必在評論區(qū)大概講一下是怎么工作的

什么是剛性系統(tǒng)

在某些條件下，使用自適應(yīng)方法求解時，會出現(xiàn)間隔不斷減少，以致需要花費(fèi)巨量算力才能得出結(jié)果，或者根本得不出結(jié)果。這種系統(tǒng)就叫做剛性系統(tǒng)。剛性與目標(biāo)方程的形狀無關(guān)，而與系統(tǒng)的特征值有關(guān)，特征值的實(shí)部越靠近負(fù)無窮剛性就越大。通常來說特征值的實(shí)部小于-2就為剛性系統(tǒng)。如何計(jì)算系統(tǒng)特征值？問得好，我也不知道

??：dy/dt=-15y就是一個典型的剛性系統(tǒng)，當(dāng)t=0時，y=1

分別使用間隔為0.25和0.125的歐拉法求得的圖像：

可以看到?jīng)]有收斂在目標(biāo)函數(shù)附近，而使用間隔為0.125的rk4方法可以得到:

如此可以看到部分方法對剛性系統(tǒng)還是有效果的

針對剛性系統(tǒng)，rk黨并沒有放棄治療，并且研發(fā)出了

隱式rk方法

顯式rk方法是從上往下計(jì)算，非常簡潔方便。而隱式rk方法是每條式子都包含自身和其他式子?

盡管隱式rk法可以普遍地求解剛性和非剛性系統(tǒng)，但是相比顯式rk方法，隱式需要的算法復(fù)雜度和算力真的大太多了。需要了解如何計(jì)算隱式rk方法的可以去查閱相關(guān)文獻(xiàn)

各種隱式rk方法：

預(yù)測-評估-校正-評估方法

實(shí)際上梯形法就是最簡單的預(yù)測評估-校正-評估方法

步驟如下:? ?1.預(yù)測: 使用歐拉法預(yù)測下一個點(diǎn)? ? 2.評估: 計(jì)算下一個點(diǎn)的斜率? ? ?3.校正: 歐拉法的斜率與評估的斜率作平均得到新斜率? ??4. 評估: 計(jì)算新斜率對應(yīng)下一個點(diǎn)的斜率 (一般來說不需要最后這一步)

這就是PECE的步驟了, 當(dāng)然最后評估得到的斜率也可以繼續(xù)拿回去糾正, 這樣就變成了PECECE方法

除非在剛性非常大的系統(tǒng)內(nèi), 這種方法會隨著糾正的次數(shù)增多而逐漸接近目標(biāo)方程. 當(dāng)做了無窮部評估校正步驟后, 必定會落在目標(biāo)方程上

應(yīng)用上述方法也是非常簡單的, 如圖:?

分別使用歐拉法, 梯形法, 和rk4方法求解著名的洛倫茲方程:? ?(間隔: .001,? 起點(diǎn):[.1, .1, .1])

紅色 - 歐拉法? ? ? ?綠色 - 梯形法? ? ? ? ?藍(lán)色 - rk4方法

下一個項(xiàng)目應(yīng)該就是隨機(jī)抽取一種幸運(yùn)方法拿去模擬三體運(yùn)動了

彩蛋：一種簡單粗暴地把數(shù)據(jù)平滑化的方法

無論在做什么項(xiàng)目，比如說上述的數(shù)值積分，總會得到一堆離散的數(shù)據(jù)，有時候就是不想用線性插值，想要各個數(shù)據(jù)直接平滑過渡，那么這里有一種懶人方法：多項(xiàng)式

如果現(xiàn)在有n組數(shù)據(jù)（這里用[ti, yi]做例子），那么假設(shè)光滑的目標(biāo)函數(shù)為一元n-1次方程：

那么可以列一個有n個方程度方程組，求出各系數(shù)就得到了完全光滑的目標(biāo)函數(shù)y(t)，不過需要注意，這個函數(shù)只有在最小t和最大t之間才好看，超過了這個范圍將會急速變?yōu)闊o窮。

什么？這么大的方程組不會解？真的沒辦法呢，列一個這樣的線性方程組：

然后讓計(jì)算機(jī)解這個就行啦