微分方程數值解法|專題（2）——偏微分方程的有限差分方法

作者：由學弱猹發表于收藏時間：2019-05-11

大家好！

趁熱打鐵，在介紹完常微分方程的那一套理論後，我們繼續介紹偏微分方程的相關內容。

學過偏微分方程的都知道，偏微分方程一般有三類：橢圓方程，拋物方程和雙曲方程。它們每個方程都各有自己的一套理論，要想完全瞭解相關性質和理論是極為困難的。但是有幸在數值解中，我們可以看到這些方程都可以用一套框架來解決。這也是這篇文章要重點關注的內容。

不過在我們進行偏微分方程的話題之前，我們還會再引入一些

之前常微分方程還沒有說完的

重要的內容和方法。因為最近我還沒有更新

有限元方法

的計劃，所以這也是專題的

最近的最後一篇

文章。

提供之前的筆記：

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

我們開始本節的內容。

線性多步法的絕對穩定區

拋物方程的有限差分方法

馮諾依曼分析法

Crank-Nicolson方法

雙曲方程的有限差分方法

橢圓方程的有限差分方法

五點格式法

線性多步法的絕對穩定區

我們還是研究方程

。無論是最簡單的向前Euler法，還是之後的龍格庫塔方法，我們的導數項

都是沒有變的。我們之前也說了它們是顯式格式。這是因為我們在計算的時候，

右端的導數項的時間點都在過去

。實際情況下，右邊的導數項是可以變的，比方說向後Euler法

$w_{i+1} = w_i + h f(t_{i+1}, w_{i+1})$

你會發現，導數項的時間是

$t_{i+1}$

，所以你會發現，這個數列通項沒有辦法再那麼容易得到了，因為你在計算時間點

$t_{i+1}$

的時候，所依賴的導數值與

$t_{i+1}$

有關。這個時候在進行程式設計的時候，一般我們採用

迭代法

。也就是說，我們把左右兩端的

$w_{i+1}$

都當作未知數

，然後做類似於這樣的迭代

$x_{n+1} = w_i + hf(t_{i+1}, x_n)$

然後做出一個關於

$\{x_n\}$

的收斂數列，不過這不是我們重點關注的內容。

既然我們可以把它改成

$f(t_{i+1},y_{i+1})$

，那麼我自然可以改成之前某些時間點的線性組合。所以我們假設

，來匯出線性多步法的一般形式

$\sum_{j =0 } ^ r \alpha_j w_{n+j} = h \sum_{j=0}^r \beta_j \lambda w_{n+j}$

設

$z = h\lambda$

，我們可以得到我們的一般形式

$\sum_{j=0}^{r}\left(\alpha_{j}-z \beta_{j}\right) w_{n+j}=0$

它對應的特徵多項式為

$R(x) = \sum_{j=0}^{r}\left(\alpha_{j}-z \beta_{j}\right) x^{j}$

，所以我們只需要考慮讓它的根滿足一定條件，即可找到它們的絕對穩定區。下面我們來舉一些實際的例子，來說明一下如何操作。

Example 1： Forward Euler Scheme

$w_{i+1} = w_i + h \lambda w_i$

兩邊寫出它的特徵多項式，可以得到

，一個一元一次方程（注意這裡的

不是變數，因為這個方程只是關於

的），所以為了要求滿足根的條件，我們就要求

，這與我們上一節推的結果完全一致。

Example 2： Backward Euler Scheme

$w_{i+1} = w_i + h \lambda w_{i+1}$

兩邊化為特徵多項式，可以得到

，解

可以得到

$x = \frac1{1-z}$

，那麼你也可以看到，只要求

$|\frac1{1-z}| < 1$

就可以，這個式子對於

的時候是恆成立的，所以向後Euler其實是無條件穩定的。

Example 3： Trapezoidal Scheme

$w_{i+1}-w_i= \frac{h}{2}(\lambda w_i + \lambda w_{i+1})$

兩邊化為特徵多項式，有

$x - 1 = \frac{z}{2}(1+x)$

，解

可以得到

$x = \frac{2+z}{2-z}$

。學過復變的話，你應該知道這個式子其實是一個保形對映。這樣的話，可以得到，如果

，那麼

$\operatorname{Re}(z) < 0$

，所以這個的絕對穩定區就是負半平面。

這裡只是一些比較簡單的例子，但是實際情況下的絕對穩定區的推導其實是相當複雜的。LeVeque的書上有一些篇幅，感興趣的朋友們可以去檢視一下。

拋物方程的有限差分方法

我們研究的就是下面的熱方程

$\left\{\begin{array}{ll}{u_{t}=c u_{x x},} & {a \leqslant x \leqslant b \quad t \geqslant 0} \\ {u(x, 0)=f(x),} & {a \leqslant x \leqslant b} \\ {u(a, t)=l(t),} & {t \geqslant 0} \\ {u(b, t)=r(t),} & {t \geqslant 0}\end{array}\right.$

學過熱方程的都知道，它本質上是一個衡量杆的溫度的方程。當然了它本質上其實是一個

初邊值問題

，直觀來看是因為，它既需要初值條件，又需要邊值條件。

下面我們來看看如何差分吧。在之前我們做過對導數的離散，那麼這裡

相當於對方程針對時間變數的導數做離散，關鍵在於

$u_{xx}$

，它是一個二階導數，那麼如何去做呢？

事實上方法是完全類似的，在之前我們用向前Euler法來進行近似的時候，其實本質上是下面這個近似。

$\frac{u_{n+1} - u_n}{h} \simeq u$

（注意我們這裡因為方程的符號改成了

，所以自然統一的用

來表示了）如果將步長標記為單位長度，那麼就相當於

注意這裡的

$\nabla$

是

向後差分

符號，

不是梯度

。所以你也能猜到，如果我們要估計

，結果是完全類似的，也就是下面的過程

我們再把步長

代回去，就可以得到我們最後的差分格式

$u_{x x}(x, t) \approx \frac{1}{h^{2}}(u(x+h, t)-2 u(x, t)+u(x-h, t))$

其中你也可以看出，這個差分格式因為是對空間變數差分的，所以時間沒有變，而空間變數中

就對應著上面的

$u_{n+1}$

。

你也可以看出，從這裡開始，因為方程開始有多個變數，所以我們不能夠再使用之前的簡單的下標標記了。這裡的

的離散也是容易的，就是下面的結果。

$u_{t}(x, t) \approx \frac{1}{k}(u(x, t+k)-u(x, t))$

因為我們在空間上設定了離散步長為

，所以我們這裡在時間上就用另外一個字母

來表示。那麼你也容易看出，分別用差分格式代到兩邊，就能得到我們最終的表示式。而且使用二維泰勒展開，我們也容易得到我們的區域性截斷誤差，分別是

$h^{2} u_{x x x x}\left(c_{1}, t\right) / 12$

和

$k u_{t t}\left(x, c_{2}\right) / 2$

，其中

$c_1 \in (x-h ,x+h), c_2 \in (t, t+h)$

，潛臺詞就是說，它的區域性截斷誤差就是

。

當然了，上一節說了那麼多穩定性，我們這一節自然也不能錯過。我們定義

$w_{ij}$

就是在

處的

數值解

，那麼自然的，我們可以將它按照時間的順序排列，得到下面的一個形式

$w_{i, j+1}=w_{i j}+\frac{c k}{h^{2}}\left(w_{i+1, j}-2 w_{i j}+w_{i-1, j}\right)$

這裡，我們假設

$\sigma = \frac{ck}{h^2}$

，那麼上面的式子就是

$w_{i j}+\frac{c k}{h^{2}}\left(w_{i+1, j}-2 w_{i j}+w_{i-1, j}\right)=\sigma w_{i+1, j}+(1-2 \sigma) w_{i j}+\sigma w_{i-1, j}$

你也可以看出來，當我只考慮時間的變化後，固定時間，再考慮空間，其實我們也可以得到一個關於空間變數的數列通項公式。因為我們的變數變成了兩個，所以這樣的一個通項就不能夠輕鬆的使用特徵多項式分析了，但是還好我們還是有招，這就是後面所使用的

馮諾依曼分析法

。

馮諾依曼（Von Neumann）分析法

雖然我們在拋物方程這裡引入了它，但實際上它可以應用在各種常見的偏微分方程中。所謂的馮諾依曼分析法，其本質就是用向量，將所有的空間離散點拼起來。那麼這樣的話，這個矩陣所對應的表示式就僅僅與時間有關了。這樣自然分析就會方便很多。

回到上面這個問題，我們針對這個差分格式，就可以化成這樣的矩陣。

$\left[ \begin{array}{c}{w_{1, j+1}} \\ {\vdots} \\ {w_{m, j+1}}\end{array}\right]=\left[ \begin{array}{ccccc}{1-2 \sigma} & {\sigma} & {0} & {\cdots} & {0} \\ {\sigma} & {1-2 \sigma} & {\sigma} & {\ddots} & {\vdots} \\ {0} & {\sigma} & {1-2 \sigma} & {\ddots} & {0} \\ {\vdots} & {\ddots} & {\ddots} & {\ddots} & {\sigma} \\ {0} & {\cdots} & {0} & {\sigma} & {1-2 \sigma}\end{array}\right] \left[ \begin{array}{c}{w_{1 j}} \\ {\vdots} \\ {w_{m j}}\end{array}\right] + \sigma \left[ \begin{array}{c}{w_{0, j}} \\ {0} \\ {\vdots} \\ {0} \\ {w_{m+1, j}}\end{array}\right]$

寫的緊湊一點，就是公式

$\mathbf{w}_{j+1} = A \mathbf{w}_j + \mathbf{s}_{j}$

。那麼為了讓這個公式收斂，我們顯然需要讓

$\rho(A)<1$

，其中

$\rho(A)$

表示它的最大特徵值。

這裡的

你也可以看到，是一個

三對角矩陣

，它是數值線性代數里非常重要的一類矩陣。這裡我們不加證明的給出一個特殊矩陣的特徵值。

Proposition 1：

對於矩陣

$\boldsymbol{T}=\left[ \begin{array}{ccccc}{1} & {-1} & {0} & {\cdots} & {0} \\ {-1} & {1} & {-1} & {\ddots} & {\vdots} \\ {0} & {-1} & {1} & {\ddots} & {0} \\ {\vdots} & {\ddots} & {\ddots} & {\ddots} & {-1} \\ {0} & {\cdots} & {0} & {-1} & {1}\end{array}\right]$

，它的特徵值為

$\lambda_{j}=1-2 \cos \pi j /(m+1), j=1, \cdots, m$

事實上，這個矩陣我們可以做一個拆分，得到

$\mathbf{T} = \mathbf{I + S}$

，其中

$\mathbf{S} = \left[ \begin{array}{ccccc}{0} & {-1} & {0} & {\cdots} & {0} \\ {-1} & {0} & {-1} & {\ddots} & {\vdots} \\ {0} & {-1} & {0} & {\ddots} & {0} \\ {\vdots} & {\ddots} & {\ddots} & {\ddots} & {-1} \\ {0} & {\cdots} & {0} & {-1} & {0}\end{array}\right]$

，根據高等代數里有關特徵值的結論，容易得到

$\mathbf{S}$

的特徵值就是

$\lambda_{j}=-2 \cos (\pi j /(m+1)), j=1, \cdots, m$

，對，就是少了一個

而已。

回到我們這個例子，你也容易透過變換得到

$A = (1-2\sigma)\mathbf{I} + \sigma \mathbf{S}$

，所以它的特徵值就是

$\lambda_j = 2 \sigma(\cos ( \pi j /(m+1) )-1)+1, j=1, \cdots, m$

。

透過對餘弦的估計，你可以知道，

$\rho(A) \in (-4\sigma+1, 1)$

，這樣的話根據

$|\rho(A)| < 1$

，就能得到我們最後的結果

$\sigma < \frac12$

。注意

我們一般不取

$\sigma = \frac12$

，這是因為數值演算法往往存在誤差。

下面我們再介紹一種

向後差分方法

。我們不再給出推導細節，直接給出它的差分格式

$\frac{1}{k}\left(w_{i j}-w_{i, j-1}\right)=\frac{c}{h^{2}}\left(w_{i+1, j}-2 w_{i j}+w_{i-1, j}\right)$

其實變化就是在左邊，把它改成了一種向後差分的格式（注意這和之前的向後Euler不是一回事）。那麼也容易透過化簡得到我們的矩陣方程。

$\left[ \begin{array}{ccccc}{1+2 \sigma} & {-\sigma} & {0} & {\cdots} & {0} \\ {-\sigma} & {1+2 \sigma} & {-\sigma} & {\ddots} & {\vdots} \\ {0} & {-\sigma} & {1+2 \sigma} & {\ddots} & {0} \\ {\vdots} & {\vdots} & {\vdots} & {\vdots} & {-\sigma} \\ {0} & {\cdots} & {0} & {-\sigma} & {1+2 \sigma}\end{array}\right] \left[ \begin{array}{c}{w_{1 j}} \\ {\vdots} \\ {w_{m j}}\end{array}\right]=\left[ \begin{array}{c}{w_{1, j-1}} \\ {\vdots} \\ {w_{m j}-1}\end{array}\right]+\sigma \left[ \begin{array}{c}{w_{0 j}} \\ {0} \\ {\vdots} \\ {0} \\ {w_{m+1, j}}\end{array}\right]$

寫的緊湊一些，其實就是

$\boldsymbol{w}_{j}=A^{-1} \boldsymbol{w}_{j-1}+\boldsymbol{b}$

。那麼我們只需要找到保證

$\rho(A)>1$

即可。這是顯然的，因為用和上面一模一樣的方法，就可以得到特徵值為

$1+2 \sigma-2 \sigma \cos \frac{\pi j}{m+1} >1$

，所以這個格式是無條件穩定的。

Crank-Nicolson方法

因為之前的兩個方法它們的誤差是

，這其實會有問題就是主要的誤差其實都在時間步長上。我們需要把時間步長

取得很小，可能才能夠達到我們要的精度。因此Crank-Nicolson方法可以很好地彌補這個缺陷。

這個方法本質上就是用不同的差分方法來近似

$u_{xx}$

和

。這裡我們的

使用向後差分，和上面那一個一樣。而我們的

$u_{xx}$

使用的是混合差分。具體來說，就是用

$\frac{1}{h^{2}}\left(\frac{1}{2}\left(w_{i+1, j}-2 w_{i j}+w_{i-1, j}\right)+\frac{1}{2}\left(w_{i+1, j-1}-2 w_{i, j-1}+w_{i-1, j-1}\right)\right)$

代替

$u_{xx}$

。如何考慮它的穩定性分析呢？如果我們還是一樣，按照時間的順序去排列，可以得到下面的形式

$-\sigma w_{i-1, j}+(2+2 \sigma) w_{i j}-\sigma w_{i+1, j}=\sigma w_{i-1, j-1}+(2-2 \sigma) w_{i, j-1}+\sigma w_{i+1, j-1}$

這個方程矩陣化會稍微麻煩一些的，大體上的形式就是

$A \mathbf{w}_j = B\mathbf{w}_{j-1} + \sigma (\mathbf{s}_{j-1} + \mathbf{s}_j)$

，其中

$A=\left( \begin{array}{ccccc}{2+2 \sigma} & {-\sigma} & {0} & {\cdots} & {0} \\ {-\sigma} & {2+2 \sigma} & {-\sigma} & {\ddots} & {\vdots} \\ {0} & {-\sigma} & {2+2 \sigma} & {\ddots} & {0} \\ {\vdots} & {\ddots} & {\ddots} & {\ddots} & {-\sigma} \\ {0} & {\cdots} & {0} & {-\sigma} & {2+2 \sigma}\end{array}\right)$

，且

$B=\left( \begin{array}{ccccc}{2-2 \sigma} & {\sigma} & {0} & {\cdots} & {0} \\ {\sigma} & {2-2 \sigma} & {\sigma} & {\ddots} & {\vdots} \\ {0} & {\sigma} & {2-2 \sigma} & {\ddots} & {0} \\ {\vdots} & {\ddots} & {\ddots} & {\ddots} & {\sigma} \\ {0} & {\cdots} & {0} & {\sigma} & {2-2 \sigma}\end{array}\right)$

。那麼現在問題自然就是落到了矩陣

$A^{-1}B$

的特徵值上。

首先，我們把矩陣做一些拆解，容易看出

$A = (2+2\sigma)\mathbf{I} - \sigma \mathbf{S}$

，

$B = (2-2\sigma)\mathbf{I} + \sigma\mathbf{S}$

。那麼這樣的話，考慮設

$\lambda_j$

為

$\mathbf{S}$

的特徵值，並且對應一個特徵向量

，那麼這就容易得到

$Bv_j = (2-2\sigma)v_j + \sigma\mathbf{S}v_j = (2-2\sigma + \sigma \lambda_j)v_j$

然後呢，如果我們有

$Av_j = \lambda_j v_j$

，那麼就一定有

$A^{-1}v_j = \frac1\lambda v_j$

。所以根據這個結論，我們容易得到最後的

$A^{-1}B$

的每一個特徵值為

$\mu_j =$

$\frac{2-2\sigma + \sigma \lambda_j}{2+2\sigma - \sigma\lambda_j}, \lambda_{j}=-2 \cos \pi j /(m+1)$

容易驗證，這裡每一個數其實都是嚴格小於1的，也就是說格式是無條件穩定的。

但是如果單純說穩定性，它並沒有什麼優勢。我們透過對它做區域性截斷誤差的誤差分析，就可以看出來更多細節。我們這裡詳細的來走一遍。這裡不失一般性，我們假設

。

首先還是考慮，使用泰勒展開，首先容易得到的是

$\frac{u(x,t) - u(x, t-k)}{k} = u_t(x, t) - \frac k2u_{tt}(x,t) + \frac{k^2}{6} u_{ttt}(x,\xi_1)$

（我希望你沒有忘記，

對應的是

$w_{ij}$

的時間點）然後呢，我們再考慮兩個中心差分，有

$\frac{1}{h^2}\left(u(x+h,t)-2 u(x,t)+u(x-h,t)\right) = u_{xx}(x,t) + \frac{h^2}{24}(u_{xxxx}(\xi_2,t) + u_{xxxx}(\xi_3,t))$

當然只要改一個時間點，就可以得到

$\frac{1}{h^2}\left(u(x+h,t-k)-2 u(x,t-k)+u(x-h,t-k)\right) = u_{xx}(x,t-k) + \frac{h^2}{24}(u_{xxxx}(\xi_4,t-k) + u_{xxxx}(\xi_5,t-k))$

如果我們略去最高階的誤差項，那麼也能夠看出，右邊最最關鍵的部分就是

$\frac{u_{xx}(x,t) + u_{xx}(x,t-k) }{2}$

，而左邊的東西就是

$u_t(x, t) - \frac k2u_{tt}(x,t)$

。我們還是抓住我們的時間和空間都在

點展開，就可以得到

$\frac{u_{xx}(x,t-k)}{2} = \frac{u_{xx}(x,t) - ku_{xxt}(x,t)}2 + O(k^2)$

剩下的就是拼在一起的事情了，這個證明已經沒什麼難度了，因為

$u_{xxt}(x,t) =\frac{\partial{u_{xx}}(x,t)}{ \partial t} = \frac{\partial{u_{t}}(x,t)}{ \partial t} =u_{tt}(x,t)$

這就完成了分析，因為你可以看到，

那一項已經被消掉了，所以這個誤差是

，這已經比上面的誤差要好很多了。

事實上，如果你把左邊的時間變數的差分離散改成中心差分，就變成了Richardson格式。這是一個無條件不穩定的格式，它的具體分析就是19年丘成桐杯競賽的應用數學個人賽的第二題。據北大小夥伴說，如果做對3個題，排名據說就能到銀牌……

雙曲方程的有限差分方法

這個時候，我們考慮的方程為

$\left\{\begin{array}{ll}u_{tt} = c^2u_{xx} \\{u(x, 0)=f(x),} & {a \leqslant x \leqslant b} \\ {u_{t}(x, 0)=g(x),} & {a \leqslant x \leqslant b} \\ {u(a, t)=l(t),} & {t \geqslant 0} \\ {u(b, t)=r(t),} & {t \geqslant 0}\end{array}\right.$

怎麼差分，我覺得你已經有數了。兩個變數都做中心差分即可，也就是格式寫成

$\frac{w_{i, j+1}-2 w_{i j}+w_{i, j-1}}{k^{2}}-c^{2} \frac{w_{i-1, j}-2 w_{i j}+w_{i+1, j}}{h^{2}}=0$

按照時間順序排一下，可以得到方程的形式為

$w_{i, j+1}=\left(2-2 \sigma^{2}\right) w_{i j}+\sigma^{2} w_{i-1, j}+\sigma^{2} w_{i+1, j}-w_{i, j-1}$

這裡的

$\sigma = \frac{ck}{h}$

。

當你把這個方程寫出來之後你可能就發現情況有點不太對了。首先就是初值的情況要單獨拉出來，其次就是這個方程即使是按照時間排序，也會與之前的

兩步時間量

有關。那應該怎麼辦呢？

我們一個一個來解決，首先初值的情況，就是這裡的

的情況，你知道

對應的就是初值的時候，這個值已經由初值條件控制住了。但是

你會發現它會與一個

的量有關係，這是個什麼玩意？也是因為這個，我們需要給出一個估計，一個關於

$w_{i, -1}$

的估計。比較好辦的是我們有了關於

的導數的表示式，所以我們可以考慮構造一箇中心差分格式來近似導數。具體點來說，我們的估計是

$w_{i,-1} \approx w_{i 1}-2 k g\left(x_{i}\right)$

代入，合併就可以得到我們的初值方程

$w_{i 1}=\left(1-\sigma^{2}\right) w_{i 0}+k g\left(x_{i}\right)+\frac{\sigma^{2}}{2}\left(w_{i-1,0}+w_{i+1,0}\right)$

事實上，它也可以被寫成是一個矩陣形式。而之後的矩陣方程都按照那個格式來寫即可，最後我們可以得到下面的矩陣方程組。

$\left( \begin{array}{c}{w_{11}} \\ {\vdots} \\ {w_{m 1}}\end{array}\right)=\frac{1}{2} \boldsymbol{A} \left( \begin{array}{c}{f\left(x_{1}\right)} \\ {\vdots} \\ {f\left(x_{m}\right)}\end{array}\right)+k \left( \begin{array}{c}{g\left(x_{1}\right)} \\ {\vdots} \\ {g\left(x_{m}\right)}\end{array}\right)+\frac{1}{2} \sigma^{2} \left( \begin{array}{c}{l\left(t_{0}\right)} \\ {0} \\ {\vdots} \\ {0} \\ {r\left(t_{0}\right)}\end{array}\right)$

$\left( \begin{array}{c}{w_{1, j+1}} \\ {\vdots} \\ {w_{m, j+1}}\end{array}\right)=\boldsymbol{A} \left( \begin{array}{c}{w_{1 j}} \\ {\vdots} \\ {w_{m j}}\end{array}\right)-\left( \begin{array}{c}{w_{1, j-1}} \\ {\vdots} \\ {w_{m, j-1}}\end{array}\right)+\sigma^{2} \left( \begin{array}{c}{l\left(t_{j}\right)} \\ {0} \\ {\vdots} \\ {0} \\ {r\left(t_{j}\right)}\end{array}\right)$

$\boldsymbol{A}=\left( \begin{array}{ccccc}{2-2 \sigma^{2}} & {\sigma^{2}} & {0} & {\cdots} & {0} \\ {\sigma^{2}} & {2-2 \sigma^{2}} & {\sigma^{2}} & {\ddots} & {\vdots} \\ {0} & {\sigma^{2}} & {2-2 \sigma^{2}} & {\ddots} & {0} \\ {\vdots} & {\ddots} & {\ddots} & {\ddots} & {\sigma^{2}} \\ {0} & {\cdots} & {0} & {\sigma^{2}} & {2-2 \sigma^{2}}\end{array}\right)$

如果我們把一般的這個方程拼的緊湊一些，我們就可以得到它的一般情形其實是諸如

$\mathbf{w}_{j+1} = A \mathbf{w}_j - \mathbf{w}_{j-1} + \sigma^2 \mathbf{s}_j$

的，所以核心還是在前面。這是一個多步的數列，它的解決方案類似於高次的常微分方程——

化為方程組

。具體點來說，我們應該把方程組寫成下面這個形式。

$\left( \begin{array}{c}{\mathbf{w}_{j+1}} \\ {\mathbf{w}_{j}}\end{array}\right)=\left( \begin{array}{cc}{A} & {-\mathbf{I}} \\ {\mathbf{I}} & {\mathbf{0}}\end{array}\right) \left( \begin{array}{c}{\mathbf{w}_{j}} \\ {\mathbf{w}_{j-1}}\end{array}\right)+\sigma^{2} \left( \begin{array}{c}{\mathbf{s}_{j}} \\ {\mathbf{0}}\end{array}\right)$

這樣我們只需要考察另一個矩陣

的特徵值就可以了。

用傳統的數值線性代數的思路當然沒有問題，但是這裡我們換一種分析方法。我們考慮設

$\lambda$

為

的特徵值，並且設

為對應的特徵向量對，那麼直接走矩陣運算，你能夠得到，

$\mu = \lambda + \frac1\lambda$

是

的特徵值。而

的特徵值我們用完全一樣的思路可以知道它應該是落在

$(2 - 4 \sigma^2 , 2)$

這個區間內。那麼可以想一想，如果

$2-4\sigma^2 < -2$

，那麼對應的

$\lambda$

是完全可以取到一個

的值，這就沒有辦法使得數值解穩定。所以我們讓

$\sigma^2 \le 1$

，因為這個時候，可以發現

$|\lambda| = 1$

，才有可能

$\mu$

取到我們的那個區間內（其實也不好說，因為兩邊習慣上是開區間，但是因為數值誤差，寫成閉區間倒也無妨）。雖然這個結果不是特別好，但是至少誤差不會呈指數增大，所以我們接受這個結果。

透過與上面相同的誤差分析，你可以得到，這個格式是

的誤差。

橢圓方程的有限差分方法

其實你也可以看出，上面的雙曲方程和拋物方程的穩定性分析方法都是使用馮諾依曼分析法的，可以說是極為方便的。但是橢圓方程的分析並不是這樣，它的方法有些另類。因此我們以它作為我們最後的結束的內容。

首先我們研究的方程一般長這個樣子

$\Delta u = f$

這裡的

$\Delta$

是拉普拉斯運算元，

不是差分符號

。考慮一下更特殊的情形，就是

$u_{xx} + u_{yy} = f$

那麼如何作差分呢？其實這個地方方法倒是類似的。不過呢我們需要給一些比較有趣的名字。

五點格式法

我們還是一樣考慮以

$w_{ij}$

來表示

處的數值解。這裡因為我們有一個額外的函式

，所以我們也需要給它一個記號

$f_{ij}$

，並且同樣的我們也假設

方向的步長為

，

方向的步長為

。有了這些假設，我們的差分方法其實就可以寫成這個形式

$\frac1{h^2}(w_{i-1, j} - 2 w _{i ,j} + w_{i+1, j}) + \frac1{k^2}(w_{i, j-1} - 2w_{i ,j} + w_{i, j+1}) = f_{ij}$

因為兩個方向其實都是空間變數，我們有理由相信它們的性質比較類似。基於這個考慮我們不妨假設

，那麼組合化簡可以得到

$\frac{1}{h^{2}}\left(w_{i-1, j}+w_{i+1, j}+w_{i, j-1}+w_{i, j+1}-4 w_{i j}\right)=f_{i j}$

下面這張圖一定程度上解釋了為什麼我們稱這樣的差分方法為“五點格式法”

也許你會想，這個表示式寫出來之後，我們也確確實實可以考慮用馮諾依曼分析法。但是問題在於，無論怎麼排布，我們都

沒有辦法

獲得諸如雙曲和拋物方程那樣好的結構。因此這邊提出了一個非常有趣的方案叫做

堆積

（stacking），通俗說就是把

所有的，每個點上的數值解排成一個列向量

。比方說我們每一個維度上取了

個點，那麼一共就會有

個點，向量自然也就是

維的了。另外排布的方式我們自然也有講究，具體來說就是這樣

$w=\left[ \begin{array}{c}{w^{[1]}} \\ {w^{[2]}} \\ {\vdots} \\ {w^{[m]}}\end{array}\right]$

，

$w^{[j]}=\left[ \begin{array}{c}{w_{1 j}} \\ {w_{2 j}} \\ {\vdots} \\ {w_{m j}}\end{array}\right]$

按照我們這樣的排列，容易發現我們的矩陣應該長這樣

$A=\frac{1}{h^{2}} \left[ \begin{array}{ccccc}{T} & {I} & { } & { } & { } \\ {I} & {T} & {I} & { } & { } \\ { } & {I} & {T} & {I} & { } \\ { } & { } & {\ddots} & {\ddots} & {\ddots} \\ { } & { } & { } & {I} & {T}\end{array}\right]$

$T=\left[ \begin{array}{ccccc}{-4} & {1} & { } & { } & { } \\ {1} & {-4} & {1} & { } & { } \\ { } & {1} & {-4} & {1} & {} \\ { } & { } & {\ddots} & {\ddots} & {\ddots} \\ { } & { } & { } & {1} & {-4}\end{array}\right]$

你也可以看出，這實際上是一個極為稀疏的大矩陣。但是還好，它的結構還算完整。而對應的方程其實就是

$A\mathbf{W} = \mathbf{F}$

，注意我們這裡的結果因為不再與時間有關，所以實際上它不再是一個迭代的數列，這一點和前面兩個有本質上的不同。

下面，我們來做一些穩定性分析吧。首先你也能看出來，因為沒有時間變數，我們這裡自然不能透過上面的計算譜半徑的方法來判斷是否穩定。所以這裡實際上問題就落在了我們怎麼處理這單獨的一個矩陣方程上。

最直觀的想法就是：考慮構造另外一個方程

$A\mathbf{W_2} = \mathbf{F_2}$

，只是說我們這裡的

$\mathbf{W_2}$

換成對應時間點的精確解。

這當然是可以的，不過這樣子得到的結果，後面計算的實際上就是精確解代入得到的一個式子。那麼如果我們相減兩式，得到

$A(\mathbf{W_2-W} )= \mathbf{F_2-F}$

，那麼左邊的

$\mathbf{W_2-W}$

就是我們的每一個點的

整體截斷誤差

。但是右邊就需要我們單獨計算了。

我們具體寫出來右邊式子的每一個元素，就是

$\tau_{i j}=\frac{1}{h^{2}}\left(u\left(x_{i-1}, y_{j}\right)+u\left(x_{i+1}, y_{j}\right)+u\left(x_{i}, y_{j-1}\right)+u\left(x_{i}, y_{j+1}\right)-4 u\left(x_{i}, y_{j}\right)\right)-f\left(x_{i}, y_{j}\right)$

相信你也知道該怎麼做了吧？這就是最經典的區域性截斷誤差的分析，它最後的結果是

$\tau_{i j}=\frac{1}{12} h^{2}\left(u_{x x x x}+u_{y y y y}\right)+O\left(h^{4}\right)$

你可以看到它的區域性截斷誤差是二階的。

把這些東西都放在右邊，就可以得到我們最後的誤差矩陣方程

$A\mathbf{E} = -\tau$

。容易看出，只要

$\|A^{-1}\|$

有界，那這個值

$\mathbf{E}$

，它的規模應該是和

$- \tau$

內的數相同的（想想為什麼），所以我們下面就需要檢查一個事情：是否

$\|A^{-1}\|$

有界。而一般來說根據範數的等價性，只要找到一個範數說明一下有界性即可。

我們還是考慮二範數，原因也很簡單，這樣子更容易分析。

首先這個矩陣是

$m^2 \times m^2$

維的，因此也就存在

個特徵值。我們不加證明的給出它的特徵值形式。

$\lambda_{p, q}=\frac{2}{h^{2}}((\cos (p \pi h)-1)+(\cos (q \pi h)-1))$

注意我們這裡的

相當於對應的

座標。比方說

$\lambda_{1,1}$

就相當於是

方向都在左（下）邊界的對應的特徵值。也就是對應

$w_{11}$

這個位置的特徵值。

那麼你也看出來，這個方程的分析不難的，因為我們有

$h = \frac1m$

，這樣的話可以看出，如果要證明

$||A^{-1}||$

有界，那麼就一定是找到原來的那個矩陣

中，

絕對值最小

的那個做分析。這裡我們自然取

$\lambda_{1,1}$

。容易得到

$\lambda_{1,1}=-2 \pi^{2}+O\left(h^{2}\right)$

這是一個有界的數，所以我們自然可以得到結論說這個格式是穩定的。

最後，我們簡單的提一下

九點格式法

。這個方法寫出來是長這樣的

$\begin{aligned} w_{i j}=\frac{1}{6 h^{2}} &\left[4 w_{i-1, j}+ w_{i+1, j}+4 w_{i, j-1}+4 w_{i, j+1}\right.\\ &+w_{i-1, j-1}+w_{i-1, j+1}+w_{i+1, j-1}+w_{i+1, j+1}-20 w_{i j} ] \end{aligned}$