微分方程數值解法|專題(1)——常微分方程的有限差分方法

大家好!

最近運氣可能比較好發現自己在丘成桐大學生數學競賽(簡稱「丘賽」)的應用數學與計算數學(Applied and Computational Math)這個track下拿到了第12名,也就是壓線入圍的意思。這可是把我這個半吊子嚇壞了哦……也是因為這個,學院安排了一些額外的學習任務(大霧),而這一系列的文章也是抓大放小,主要關注方法論和實用性。諸如微分方程數值解法,其實是有套路可循的,而掌握這一種套路往往就能夠用來解決很多很多問題。所以這篇專題我們將關注數值解的相關方法,並給出相關的應用舉例。

微分方程數值解法這一個學科需要的前置知識,除了數分高代以外,還需要數值分析和常微分方程,偏微分方程的一套理論,所以確實對於初學者來說已經完全不能算是「零門檻」了。所以這一系列專題,我們不會特別的關注理論(因為我不可能為了這個專門搬出一本數值分析)。相反,我相信直截了當的介紹方法,反而會讓那些不是特別鑽理論的應用者們有所受益。

我們follow的課本是Timothy, SauerNumerical Analysis(好吧,我還是搬出來了一本……),李榮華的《微分方程數值解法》,LeVequeNumerical Methods of Partial Differential Equations 以及我們學院許傳炬老師的課程材料。因此在任何單獨一本教材上,你都無法找到這份筆記上的所有內容。但是也因為是專題的緣故,參考的書本較多,儘管我在寫作前已經事先組織過,但確實可能邏輯線上不如查只查閱一本教材那麼縝密。還請大家諒解。

廢話就說到這裡,我們開始正題吧。

目錄

  • 一階常微分方程初值問題的有限差分方法與誤差分析
    • 向前Euler法及誤差分析
    • 梯形法
    • 龍格庫塔方法簡介
  • 穩定性
    • 絕對穩定區
  • 剛性常微分方程的有限差分方法
    • 龍格庫塔——切比雪夫方法

一階常微分方程初值問題的有限差分方法與誤差分析

在這裡,我們主要分析的是下面這個方程

left{egin{array}{l}{y^{prime}=f(t, y)} \ {y(a)=y_{a}} \ {t in[a, b]}end{array}
ight.

你可以看到,它的導數是一個函數 f ,我們稱它為「初值問題」的原因是,這個方程本質上描述的是一個物體在給定導數(也就是速度)後,隨著時間的變化趨勢。也就是說它主要是與時間變數有關。

所謂數值解法,本質上就是一種離散。我通過對時間步長作分割,通過一種特殊的數列迭代的公式,使得在我給定數列的前幾項之後,通過數列的迭代公式得到我每一個分點上的數值解。因為很多時候,我們需要用到方程的解在某些點上的值,而傳統的偏微分方程,很多時候都沒有辦法找出顯式解,這就是數值解所起作用的地方了。

你也可以看到,數值解必然不可能和精確解完全相同的,會產生誤差。我們下面給出與它有關的誤差定義。首先很明顯的就是整體截斷誤差,它一般定義為

Definition 1: Global Truncation Error

g_i = |w_i - y_i|

其中 y_i 就是我們的精確解。而 w_i 就是給定計算格式(也就是數列的迭代公式)之後,你所計算出來的數值解。

第二個誤差是局部截斷誤差,我們之後來解釋為什麼我們需要這個概念。它的定義是

Definition 2: Local Truncation Error

e_{i+1} = |w_{i+1} - z(t_{i+1})|

首先我們需要解釋, z(t_{i+1}) 是下面這個初值問題的解。

left{egin{array}{l}{y^{prime}=f(t, y)} \ {yleft(t_{i}
ight)=w_{i}} \ {t inleft[t_{i}, t_{i+1}
ight]}end{array}
ight.

你可以看到,這個初值的時間是從 t_i 開始的,它相當於是一小段時間內的 y 的表現。值得關注的是,它的初值設成了 w_i ,也就是說它本質上的意思就是,當我的初值精確解就是我的數值解的時候,我在下一個時間點的值應該滿足的方程。所以我們這個時候假設的誤差,其實是相當於假設之前的計算都是精確的情況下,下一步計算產生的誤差。這也是「局部」的含義。一般情況下,「局部」帶來的誤差是容易估計的,這樣的話,整體的誤差就可以通過局部的誤差來累加得到,所以引入這個局部截斷誤差,多是為了之後的誤差分析方便。

有了這些誤差的定義,我們就可以來分析相關的計算格式了。

向前Euler法及誤差分析

我們先介紹向前Euler方法。在此之前,我們設時間是等步長的,並且步長為 h

Definition 3: Forward Euler Scheme

egin{array}{l}{w_{0}=y_{0}} \ {w_{i+1}=w_{i}+h fleft(t_{i}, w_{i}
ight)}end{array}

這樣子直接看可能會糊,我們其實稍稍修改一下,可以得到下面的這個結果。

f(t_i,w_i) = frac{w_{i+1}-w_i}{h}

如果你把這些數值解都用精確解去代替的話,式子就會變成這樣

f(t_i,w_i)simeq frac{y(t_{i+1},w_i)-y(t_i,w_i)}{h}

你可以看到,當我令 h 	o 0 的時候,右邊其實就是導數的定義表達式。所以數值格式的本質,就是使用我的解的線性組合,去估計導數。而如何估計導數其實方法各有不同,而這也自然產生了不同的格式。

它的誤差好分析嗎?其實都是一個套路,這個套路叫做泰勒展開。我們用它舉例。

注意到 t_{i+1} = t_i + h ,所以通過泰勒展開,我們可以得到

yleft(t_{i}+h
ight)=yleft(t_{i}
ight)+h y^{prime}left(t_{i}
ight)+frac{h^{2}}{2} y^{prime prime}(c)

其中 c in (t_i,t_{i+1}) 。你可以看到,這裡的 y(t_i)w_i 是不一定相等的。所以如果直接推導總體的誤差,其實是會有些麻煩的(不過麻煩的點並不是這兩個,而是之後的 y(t_i)f(t_i,w_i) )。所以我們退而求其次,先求局部截斷誤差,再考慮之後的分析。

我們在前面提到說,局部截斷誤差的本質就是假設之前的數值解都是精確的。因此這裡相當於 y(t_i) =w_i 。這樣的話,我們可以得到

yleft(t_{i+1}
ight)=w_{i}+h fleft(t_{i}, w_{i}
ight)+frac{h^{2}}{2} y^{prime prime}(c)

這個時候,這裡的 y(t_{i+1}) 其實就對應了上面的 z(t_{i+1})

我們自然還可以把數值格式寫出來,得到下面的式子

w_{i+1}=w_{i}+h fleft(t_{i}, w_{i}
ight)

我們把兩式相減就可以得到 e_{i+1}=left|w_{i+1}-yleft(t_{i+1}
ight)
ight|=frac{h^{2}}{2}left|y^{prime prime}(c)
ight| le frac{Mh^2}{2}

我們有了局部的截斷誤差,自然需要考慮整體截斷誤差。剛開始顯然是沒問題的,然後我們做了一步之後,得到的整體截斷誤差其實就和局部截斷誤差沒差,關鍵是第二步之後。因為第二步之後,除了局部截斷誤差,你在第一步得到的整體截斷誤差也會變大,因此這實際上是好幾個因素的綜合。所以我們需要拆分來看。

我們定義 z(t_2) 是假設 t_1 時間點精確的情況下得到的方程的精確解。那麼顯然如果我們需要估計 g_2 = |w_2 - y _2| ,最直接的方式就是補一個 z(t_2) 進去。因為 |w_2 - z(t_2)| 我們是知道的,就是局部截斷誤差。所以不知道的是 |z(t_2) - y_2| 。我們先把這些推導寫在下面。

g_{2}=left|w_{2}-y_{2}
ight|=left|w_{2}-zleft(t_{2}
ight)+zleft(t_{2}
ight)-y_{2}
ight|leqslantleft|w_{2}-zleft(t_{2}
ight)
ight|+left|zleft(t_{2}
ight)-y_{2}
ight|

為了估計出這個式子,我們回想一下常微分方程的一些要求。在考慮常微分方程的連續依賴性的時候,我們有給出過方程所需要滿足的Lipschitz條件。我們寫在下面。

Definition 3: Lipschitz Condition

設存在 L 為常數,使得 f(t,y)S = [a, b] 	imes [y_1, y_2] 中任意的 (t, y_1), (t, y_2) 都有 |f(t,y_1) - f(t, y_2) | le L |y_1 - y_2| ,則稱 f(t,y)S 中關於變數 y 滿足Lipschitz條件。

我們引入這個條件的原因相信你也能看出來:通過導數的差距,衡量出在不同的初值下,方程行為的差異程度。其實如果你的眼睛夠尖,你就能看出來,雖然 z(t_2)y_2 多了一個假設,但是它們的微分方程的形式是完全相同的,只是初值不同。有了這個思路,我們給出下面的性質

Proposition 1:

假設 Y(t),Z(t) 是對應微分方程 y^{prime}=f(t, y) 的不同初值,那麼有 |Y(t)-Z(t)| leqslant mathrm{e}^{L(t-a)}|Y(a)-Z(a)|

簡單說明一下這個結論。設 u(t) = Y(t)-Z(t) ,那麼它顯然要不是一直正,要不是一直負的。這樣的話,我們可以去掉兩邊的絕對值。再考慮 u(t) ,這樣的話就會有下面的結果

|u^{prime}|=|f(t, Y)-f(t, Z)| leqslant L|Y(t)-Z(t)|=L|u(t)|=L u(t)

也就是說 u le  Lu(t) (左邊去掉絕對值就好),再根據中值定理可以得到

frac{ln u(t)-ln u(a)}{t-a} leqslant L

化簡即可得到我們最後的結果。

使用這個估計,我們再回頭來看 |z(t_2) - y_2| 。你可以看到,其實 z(t_2) 就是比 y_2 變了個初值而已。那麼為什麼這個定理有用呢?你看上面的定理,其實就相當於說,如果我們把時間從 t 「退回」到 a ,那麼就會使得誤差降低一個膨脹因子 e^{L(t-a)} 。所以我們有下面的結果

left|w_{2}-zleft(t_{2}
ight)
ight|+left|zleft(t_{2}
ight)-y_{2}
ight| le  e_{2}+mathrm{e}^{L h} g_{1}=e_{2}+mathrm{e}^{L h} mathrm{e}_{1}

首先,我們考慮的是第二步,那麼往回退一步(這裡也就是退了一個步長 h ), z(t_2) 就退回到了 z(t_1) ,而 y_2 就退回到了 y_1 。所以你可以看到,這個時候兩個值的差就恰好是第一步的整體截斷誤差。而第一步的整體截斷誤差恰好就等於第一步的局部截斷誤差。所以當我們假設了之前的數值解都精確的時候,我們其實只能夠給出第一個式子 |w_2 - z(t_2)| 的誤差估計,也就是局部截斷誤差,而其它的誤差(這裡指的是 left|zleft(t_{2}
ight)-y_{2}
ight| )其實是從之前步的整體截斷誤差再擴大所得到的。

有了這個思路,相信你也不難做後面的推導,我們直接寫出一般的結果

g_{i}=left|w_{i}-y_{i}
ight| leqslant e_{i}+mathrm{e}^{L h} e_{i-1}+mathrm{e}^{2 L h} e_{i-2}+cdots+mathrm{e}^{(i-1) L h} e_{1}

所以通過這個方法,我們就將整體的誤差,用一系列局部截斷誤差估計出來了。

如果我們進一步假設 e_i le Ch^{k+1} (我們之後說假設這個幹啥),就可以得到下面的一些不等式估計。我就直接抄書了……

g_{i} leqslant C h^{k+1}left(1+mathrm{e}^{L h}+cdots+mathrm{e}^{(i-1) L h}
ight)=C h^{k+1} frac{mathrm{e}^{i L h}-1}{mathrm{e}^{L h}-1}

leqslant C h^{k+1} frac{mathrm{e}^{Lleft(t_{i}-a
ight)}-1}{L h}=frac{C h^{k}}{L}left(mathrm{e}^{Lleft(t_{i}-a
ight)}-1
ight) (注意 ih le t_i -a

所以你可以看出來,對於這種方程,如果局部截斷誤差是 k+1 階的,就可以得到整體截斷誤差是 k 階的。

梯形法

考慮到我自己在開始分析梯形法的格式時也是處處碰壁,這裡還是決定提一下這個格式。

Definition 4: Trapezoidal Scheme

egin{array}{l}{w_{0}=y_{0}} \ {w_{i+1}=w_{i}+frac{h}{2}left(fleft(t_{i}, w_{i}
ight)+fleft(t_{i}+h, w_{i}+h fleft(t_{i}, w_{i}
ight)
ight)
ight)}end{array}

你看著感覺挺複雜的,但其實這第二個方程只是取了兩個時間的中值而已。

為了筆記的多樣性,我們用另外一個方法來解釋這個格式的來源。其實你只需要移個項,得到

w_{i+1}-w_{i} =frac{h}{2}left[fleft(t_{i}
ight)+fleft(t_{i}+h
ight)
ight]

(為了方便我們省去了第二個參數)然後我們還是,用精確解去代替這裡的數值格式,得到

y(t_{i+1})-y(t_{i}) simeq frac{h}{2}left[fleft(t_{i}
ight)+fleft(t_{i}+h
ight)
ight]

這個估計為什麼合理呢?除了用導數的估計來看,還可以通過數值積分來看。首先注意到 y(t_{i+1} ) - y(t_i) = int_{t_{i}}^{t_{i}+h} y(t) mathrm{d} t = int_{t_{i}}^{t_{i}+h} f(t) mathrm{d} t ,那麼這個時候,你可以看到,右邊的式子其實就是一個數值積分格式(梯形格式),因此這個格式也被稱為梯形格式。

我們再走一遍流程,做一遍誤差分析。你可以看到上面其實我們就是重點關注了兩個式子:精確解的泰勒展開,數值解的迭代格式。所以我們這裡也是一樣的思路,首先精確解的泰勒展開我們不難得到,可以寫成

y_{i+1}=y_{i}+h y^{prime}left(t_{i}
ight)+frac{h^{2}}{2} y^{prime prime}left(t_{i}
ight)+frac{h^{3}}{6} y^{prime prime prime}(c)

當你寫到這個式子的時候其實你就應該會感覺分析上會有些麻煩了。數值格式上既有 f(t_i) 又有 f(t_i+h) ,展開的話究竟在哪一個點展開?又怎麼樣與精確解的泰勒展開式子相消?所以我們這裡最關鍵的,自然就是對 y 做處理,這就需要我們下一個式子了。

y^{prime prime}(t)=frac{partial f}{partial t}(t, y)+frac{partial f}{partial y}(t, y) y^{prime}(t)=frac{partial f}{partial t}(t, y)+frac{partial f}{partial y}(t, y) f(t, y)

有些複雜,但它就是鏈式法則而已。

既然你發現了我們使用的本質上是二維的泰勒展開,那麼相同的思路自然也可以用到我們的數值格式上。顯然我們推導局部誤差是方便的,所以只需要假設 w_i = y_i ,那麼我們會得到下面的結果

w_{i+1}=y_{i}+frac{h}{2}left(fleft(t_{i}, y_{i}
ight)+fleft(t_{i}+h, y_{i}+h fleft(t_{i}, y_{i}
ight)
ight)
ight)

顯然第二個式子可以用相同的方法作展開,也就是說

fleft(t_{i}+h, y_{i}+h fleft(t_{i}, y_{i}
ight)
ight) = fleft(t_{i}, y_{i}
ight)+hleft(frac{partial f}{partial t}left(t_{i}, y_{i}
ight)+fleft(t_{i}, y_{i}
ight) frac{partial f}{partial y}left(t_{i}, y_{i}
ight)
ight)+Oleft(h^{2}
ight)

這樣的話,其實可以看到,再兩項作相減之後,都只會有 3 階項存在,所以可以認為局部截斷誤差 e_{i+1} = O(h^3) 。那麼整體方法就是一個 2 階方法。

關於隱式格式,我們不再做詳細的介紹,雖然我們之後還是會用到它們。

龍格庫塔(Runge-Kutta)方法簡介

這個方法本質上是用導數在各個點的取值來構造數值積分,估計函數的端點值之差。也就是估計

y(t+h) - y(t) = int_t^{t+h}f(	au, y(	au))d	au

右邊的這個積分。龍格庫塔法的意思就是:取定 t_1 = t le t_2 le t_3 le cdots le t_m le t+h ,設 k_i = f(t_i, u(t_i)) ,那麼想法就是找到一組係數 {c_i} ,使得 sum_{i=1}^{m} c_ik_i simeq f ,也就是儘可能的去逼近導數。

一個顯然的問題是:如何給出 k_i ?這個確實是沒有辦法的,因為它可能會和 y(t_i) 有關,而這個你一般是不知道的。那麼這個時候,我們給出了這麼一個逼近:先給出

(t_1, k_1) = (t_1, f(t_1,y(t_1)))

(左邊界上的點的值我們知道),然後對於 t_2 時間,用向前Euler估計,得到下面的估計

k_2 simeq f(t_2, y(t_1) + (t_2-t_1)k_1)

k_3 simeq f(t_3, y(t_1) + (t_2 - t_1)k_1 + (t_3-t_2)k_2)

(注意他們一步一步都是約等於過來的)問題就落在了:如何取這些係數?

那麼我們推廣上面的思想,可以得到下面的這個推導公式。

egin{array}[l]  kk_1 = f(t, y) \  k_2 = f(t+ha_2, y(t) + hb_{21} k_1), b_{21} = a_2 \ k_3 = f(t+ha_3, y(t) + h(b_{31}k_1 + b_{32} k_2)), b_{31} + b_{32} = a_3 \ cdots \ k_m = f(t + ha_m, y(t) + h sum_{j=1}^{m-1} b_{mj}k_j), sum_{j = 1}^{m-1} b_{mj} = a_m end{array}

其中 t_i = t + a_ih

為什麼說這是我們上面的那個思路的一種推廣呢?這是因為我們這裡其實每一步都運用了前面的信息,並且要估計未知的 y(t_i) 以便我們計算 k_i 。那麼你也可以看到,這個時候呢,我們自然需要的是:當 f(t,y) 中的時間參數 t 在某一個點 t_i 時,它的函數值參數 y 也必須要是 y(t_i) 的一個線性估計。這就是 sum_{j = 1}^{n-1} b_{nj} = a_n 的含義(不然時間就不同步了)。

然而,這只是最最一般的情況,我們很顯然不能就到此為止了。因為我們要逼近導數,所以我們考慮把真值的表達式寫成下面的這種形式

y(t+h) = y(t) + hvarphi_T(t, y, h)

而這裡的 varphi_T 是一個表達式,因為這裡的 y 都是精確解,所以這個表達式是可以通過泰勒展開得到的。很明顯,你想要幾階的精度,就把它展開到幾階即可。但是要注意的是每多展開一階,計算的複雜度就會多一個量級。比方說我們把它展開到三階,那麼表達式就可以寫成下面這樣

y(t+h) = y(t) + sum_{l = 1}^3 frac{h^l}{l!}y^{(l)}(t) + O(h^4) = y(t) + h varphi_T(t, y, h)

其中

egin{cases}  varphi_T (t, y, h) = f + frac12 h F + frac16 h^2 (F f_y + G) + O(h^3) \ F  = f_t + ff_y \  G = f_{tt} + 2 f f_{ty} + f^2 f_{yy} end{cases}

你就會發現其實這個係數就已經相當複雜了。

那麼精確解是這樣的一個形式,對應的數值格式呢?其實就是把公式改成了下面這個形式

w_{i+1} = w_i + h varphi(t, y, h)

這裡的 w_i 就是時間為 t 時刻的數值解, w_{i+1} 啥意思自然不必多說。你可以看到其實本質上的差異就是 varphivarphi_T 。而這個 varphi 其實就是我們之前假設的那個 sum_{i=1}^{m}c_ik_i, sum_{i=1}^{m} c_i = 1 。這個計算事實上也不是特別容易的,如果你要求 m 的數量比較大的話。我們這裡隨便計算兩個表達式如下

k_1 = f(t, u) = f

egin{array} kk_2 &= f(t + ha_2, u + h a_2 k _1 ) \ & = f + h a_2 F + frac12 h^2 a_2^2 G + O (h^3) end{array}

其中 F,G 參考上面的表達式。所以比方說你要求 m=2 ,那你就需要把 k_1,k_2 計算出來,但是你也能看出來這個計算量已經是非常大了。

最後,我們所需要做的事情,就是在確定了我們要的誤差階數之後,根據不同的 m ,去匹配不同的係數,使得 varphi_T, varphi 在係數上匹配。這便是龍格庫塔方法的全部思想,而它的具體格式實現,大家可以通過wikipedia尋找更多內容。

最後,提一下,我們對於龍格庫塔方法其實是要控制兩個參數 m,p ,而這對應得到的格式,我們就稱他們為 mp 階龍格庫塔方法。

穩定性

我們用剛性的常微分方程來看我們為什麼需要穩定性。一個比較簡單的方程形式是

y = epsilon y

其中 epsilon 是一個絕對值很大的負數(比如說 -10^{5} )。你自然會疑惑,這也算是上面的方程的一個例子,我們如果使用向前Euler方法,你可以知道它具有一階的精度(誤差為 O(h) ),也就是說你肯定不怕它不收斂。既然收斂了,為什麼還要擔心其它的東西?

這就牽涉到一個概念叫做穩定性。也就是說,這個格式確實收斂沒有錯,但是很有可能因為你的其它參數(比如步長)選的不好,導致它在給定的步數下沒有辦法計算到給定的精度。所以我們先介紹穩定性的相關概念,再在之後觀察為什麼剛性問題會有其它需要注意的地方。

絕對穩定區

我們還是分析之前分析過的常微分方程,但是稍微寫的具體一點,像這樣

y^{prime}(t)=lambda y(t)+g(t)

如果我們使用向前Euler法進行計算,我們就可以得到下面的計算格式

w_{i+1} = w_i + h (lambda w_i + g(t_i)) = (1+hlambda)w_i + h g(t_i)

同樣的,我們對精確解的式子做泰勒展開,可以得到

y_{i+1} = y_i + h y_i + frac{h^2}{2} y(c) = y_i + h(lambda y_i + g(t_i)) + 	au_n

其中 	au_n 就是式子最後的那個誤差項。相減之後就可以得到我們下面的這個誤差式子

E_{i+1} = (1+hlambda)E_i - 	au_n

其中這裡的 E_i = w_i - y_i (我們沒有再加絕對值)。所以你可以看出來,儘管我們存在一個局部的誤差,但是如果這裡的膨脹因子 1+hlambda 過大,你就自然會發現,一步一步的迭代,其實誤差是會發散的,並不是我們之前想的「收斂」。這可能會與之前的推導矛盾,你會覺得,明明我之前都找到階數了,結果反而發散了,這不是胡扯嗎?但是要注意的是,我們給出階數,其實只是相當於給出了誤差與步長的關係。只有你的步長趨於0的時候,我們才能得到收斂性。那麼如果你的 1+hlambda 過大,一般情況下都是 h 取得太大了,那自然不存在收斂這一說。

運用幾乎相同的分析方法,你可以觀察一下,如果回到上面一般的方程 y = f(t,y) ,那麼Euler格式一般就是有諸如

w_{i+1} = w_i + k f(w_i)

的形式,這個形式的分析方法和上面類似的,但是這裡的膨脹因子變成了 1+hL ,這裡的 L 就是函數 f 的Lipschitz常數,也就是取使得

|f(x)-f(y)| le L|x-y|

的最小的常數 L

那麼根據這些特殊的例子,你也能看出來如何定義絕對穩定區了。

Definition 5: Absolute Stability Area

定義特徵多項式 f(z), z = hlambda ,那麼滿足 |f(z)| < 1 的區域稱為絕對穩定區。

這個定義是很一般的,不同的數值格式會對應有不同的特徵多項式。比如說上面的向前Euler法,它的特徵多項式就是 1+z 。也就是 y_i 之前那個我們說的「膨脹因子」。所以你也能看出來,這個格式的絕對穩定區其實就是 |1+z| < 1 ,它在複平面上就是一個圓,當然了如果限制在實數上,對應的 z 的取值就是 (-2,0)

在這部分最後我們簡單提一下數值格式的各種可能的穩定性。

Definition 6: Stability

若對於任意的步長 h ,數值格式計算得到的解都會收斂,則稱為無條件穩定。若需要步長滿足一定的條件才可穩定,則稱為條件穩定。否則稱為無條件不穩定。

一些具體的例子會牽涉到多步法的概念,但是多步法的解決方法和框架等其實應該是偏微分方程數值解里更為常見的,所以我們這裡不會介紹這樣的方法,但是我們還是用一個具體的多步法例子,來介紹一般的線性差分方程,如何判斷它的穩定性

Example 1:

w_{i+2} - 3w_{i+1} + 2w_i = -h f(w_i)

這是一種多步法,在進行之前,我們先看看這個格式為什麼是存在的。老方法,用精確解去替代數值解,我們可以得到

y(t_{i+2}) - 3y(t_{i+1}) + 2y(t_i) simeq -hf(t_i, y(t_i))

也就是說

frac{-y(t_{i+2}) + 3y(t_{i+1}) - 2y(t_i)}{h} simeq f(t_i, y(t_i))

我們對兩個不是在時間 t_i 的式子做一下泰勒展開,看一下我們的結果

y(t_{i+2}) = y(t_i) + 2hy(t_i) + 2h^2 y(c_1)

y(t_{i+1}) = y(t_i) + h y(t_i) + frac{h^2}{2}y(c_2)

所以代入之後,你會發現 y(t_i) 的項都被消掉了,而 y(t_i) 項之前的係數就是1,和右邊的正好匹配上,剩下的就是一個小量了。所以這個式子是一個合理的導數的估計。

好的,那我們如何判斷呢?這裡給大家介紹一個數列的特徵方程法。如果你是競賽黨出身,你絕對不可能對這個名詞感到陌生。

針對這個數值例子,我們忽略所有的與右端導數項有關的內容,只考察差分的數列項。那麼針對這個例子,它對應的特徵方程其實就是 x^2-3x+2=0 ,可以得到它的兩個根為 x = 1, 2 。這樣的話,理論證明了它的根的形式就是 w_i = c_1 + c_22^i 。你看這個 2^i 就可以說明了這個格式其實是無條件不穩定的,我們不應該使用它。

一般的,我們有下面的性質:

Proposition 2:

對於特徵方程,如果它的根絕對值均小於1,則格式無條件穩定。若根均絕對值大於1,或者存在一個以上的根絕對值為1,則格式無條件不穩定。如果存在一個根絕對值為1,其餘根絕對值小於1,則格式條件穩定。

具體為什麼,我們在之後再說。

剛性常微分方程的有限差分方法

有了穩定性的概念,我們再回頭看方程 y = epsilon y 。如果我們用向前Euler法做離散,可以得到它的對應的值 z = h epsilon ,那麼你容易發現,若要求 |1+hepsilon| < 1 ,其實對於 h 的要求是相當高的(如果 epsilon = 10^{-5} ,那麼 h in (0, 2 	imes 10^{-5}) )。這個時候迭代的話,在機器上受到誤差的影響,其實很多時候會嚴重影響數值解的有效性。究竟如何修改格式和解,其實是一個大的主題。我們這裡只介紹部分方法和背後的思想,作為這一節的最後一部分。

我們之前提到了絕對穩定區的概念。那麼你也能看到,如果它的絕對穩定區涵蓋整個複平面的左半部分,那麼無論我的 h 怎麼取,在這個例子中它都是穩定的。這個格式是存在的,比如說向後Euler格式,它的穩定性要求是 frac{1}{|1-z|} < 1 ,所以至少 z<0 的時候它絕對是可以滿足的。而 z = hlambda < 0 其實是一個非常松的要求,因為 lambda <0 (如果 lambda>0 ,那麼你也能看到,這個方程本身就沒有辦法得到一個收斂的數值解,這不是數值解法喜歡關注的內容)。

所以你也可以看到,如果它的穩定區範圍很廣,那麼自然剛性常微分方程就能解。但是這種方法多半是隱式方法,而隱式方法因為通過迭代法解初值,所以計算複雜度要比顯式方法大得多。所以有沒有顯式的格式,能夠用於解剛性方程?答案是肯定的。這就是我們下面要介紹的內容。

龍格庫塔——切比雪夫方法

首先我們回憶一下龍格庫塔方法(上面雖然寫了,但是可能比較抽象?)。一個具體的例子如下:

Example 2:

egin{array} yy^* = y_i + frac12 h f(y_i, t_i) \ y_{i+1} = y_i + h f(y^*, t_i + frac12 k) end{array}

我們可以看出它就是取了一個 t_i,t_{i+1} 的中點。然後我直接用了中間點的導數值去近似估計了 t_i 點的導數值。

如果我們設 y = lambda y ,代入我們可以得到, y^* = (1 + frac12 lambda h) y ,而

y_{i+1} = y_i + h lambda y^* = y_i + hlambda ( 1 + frac12 h lambda) y_i = (1 + hlambda + frac12 (hlambda)^2)y_i

因此歸納法容易得到,如果我們的方法是 r 步的,就可以將我們的數值格式寫成這樣的形式

y_{i+1} = R(z) y_i

其中 z = hlambda ,而 R(z) = sum_{i=0}^{r}d_iz^i

如果我們需要它具有一階精度,這裡我們顯然需要 d_0 = 1, d_1 = 1 。那麼你可以看到,這個要求等價於要求多項式滿足 R(0) = 1, R(0) = 1 。學過數值逼近的同學肯定已經開始有感覺了。這可以通過正交多項式來擬合。因此我們的這個方法帶上了切比雪夫的名字。

我們不加證明的給出我們考慮的多項式

R(z) = T_r(1 + z/r^2)

容易發現,如果要求格式穩定,那麼必須要有 |R(z)|<1 。而切比雪夫多項式的性質告訴我們,它的滿足的穩定區間範圍為 (-2r^2, 0) ,你也可以看到,通過改變 r ,那麼即使是顯式的方法,我們也能夠達到我們想要的穩定性結果。

當然,如果存在擾動,那就又會有很多其它的問題,書上提到了一個位移切比雪夫多項式法 (Shifted Chebyshev Polynomial),礙於篇幅我們不再介紹,感興趣的可以尋找相關資料,或直接參考LeVeque的書。

小結

和《筆記整理》系列文章不同,《專題》系列文章會顯得有點「左一塊右一塊」的意思。而且因為我們更加註重於直接性,用於鋪墊的語句也做了大量的刪減。但還好,我們算是用了不短的篇幅,介紹完了大部分常微分方程的有限差分方法。包括顯式格式,龍格庫塔方法,穩定性分析和擴展的剛性常微分方程的解決方案。

關於隱式格式和與穩定性有關的更多概念,雖然在常微分方程數值解法中就已經出現,但是它們的應用更多的其實是在偏微分方程的差分方法中。因此我們會在之後再提到這些概念。

——————————————————————————————————————

本專欄為我的個人專欄,也是我學習筆記的主要生產地。任何筆記都具有著作權,不可隨意轉載和剽竊

個人微信公眾號:cha-diary,你可以通過它來有效的快速的獲得最新文章更新的通知。

專欄目錄:筆記專欄|目錄

想要更多方面的知識分享嗎?歡迎關注專欄:一個大學生的日常筆記。我鼓勵和我相似的同志們投稿於此,增加專欄的多元性,讓更多相似的求知者受益~


推薦閱讀:

TAG:微分方程數值解 | 微分方程 | 計算數學 |