馮諾依曼穩定性分析

數值分析中, 馮諾依曼穩定性分析 (亦作傅立葉穩定性分析) 用於驗證計算線性偏微分方程時使用特定有限差分法的數值穩定性^[1]，該分析方法基於對數值誤差的傅立葉分解。1947年英國研究人員約翰·克蘭克和菲利斯·尼科爾森在文章中對該方法進行了簡要介紹^[2]， 爾後又出現在馮諾依曼合作的文章中 ^[3] 。 洛斯阿拉莫斯國家實驗室對該方法進行了進一步發展。

數值穩定性

數值穩定性與數值誤差密切相關。使用有限差分方法進行計算時，若任意時間步的誤差不會導致其後計算結果的發散，則可稱該有限差分法是數值穩定的。如果誤差隨著進一步計算降低最終消失，該算法被認為穩定；若誤差在進計算中保持為常量，則認為該算法「中性穩定」。但如果誤差隨著進一步計算增長，結果發散，則數值方法不穩定。數值方法的穩定性可以通過馮諾依曼穩定性分析得到驗證。穩定性一般不易分析，特別是針對非線性偏微分方程。

馮諾依曼穩定性方法只適用於滿足 Lax–Richtmyer 條件 （Lax 等價定理） 的某些特殊差分法： 偏微分方程系統須線性，常係數，滿足周期性邊界條件，只有兩個獨立變量，差分法中最多含兩層時間步^[4]。 由於相對簡單，人們常使用馮諾依曼穩定性分析代替其他更為詳細的穩定性分析，用以估計差分方法中對容許步長的限制。

方法描述

馮諾依曼誤差分析將誤差分解為傅立葉級數。為了描述此過程，考慮一維熱傳導方程

${\displaystyle {\frac {\partial u}{\partial t}}=\alpha {\frac {\partial ^{2}u}{\partial x^{2}}}}$

空間網格間隔為 ${\displaystyle L}$, 對網格作 FTCS （Forward-Time Central-Space，時間步前向歐拉法，空間步三點中心差分) 離散處理，

${\displaystyle \quad (1)\qquad u_{j}^{n+1}=u_{j}^{n}+r\left(u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}\right)}$

其中 ${\displaystyle r={\frac {\alpha \,\Delta t}{\Delta x^{2}}}}$ 。${\displaystyle u_{j}^{n}}$ 為離散網格上的數值解，用於近似此偏微分方程的精確解 ${\displaystyle u(x,t)}$ 。

定義捨入誤差 ${\displaystyle \epsilon _{j}^{n}=N_{j}^{n}-u_{j}^{n}}$ 。 其中 ${\displaystyle u_{j}^{n}}$ 是離散方程 (1) 式的精確解，${\displaystyle N_{j}^{n}}$ 為包含有限浮點精度的數值解。 因為精確解 ${\displaystyle u_{j}^{n}}$ 滿足離散方程, 誤差 ${\displaystyle \epsilon _{j}^{n}}$ 亦滿足離散方程 ^[5]：

${\displaystyle \quad (2)\qquad \epsilon _{j}^{n+1}=\epsilon _{j}^{n}+r\left(\epsilon _{j+1}^{n}-2\epsilon _{j}^{n}+\epsilon _{j-1}^{n}\right)}$

此式將確定誤差的遞推關係。方程 (1) 和 (2) 中，誤差和數值解隨時間具有一致的變化趨勢。對於含周期性邊界條件的線性微分方程，間隔 ${\displaystyle L}$ 上的空間部分誤差可展開為傅立葉級數

${\displaystyle \quad (3)\qquad \epsilon (x)=\sum _{m=1}^{M}A_{m}e^{ik_{m}x}}$

其中波數 ${\displaystyle k_{m}={\frac {\pi m}{L}}}$，${\displaystyle m=1,2,\ldots ,M}$， ${\displaystyle M=L/\Delta x}$。 通過假設誤差幅度 ${\displaystyle A_{m}}$ 是時間的函數，可以給出誤差和時間的關係。 不難知單步中，誤差隨時間指數增長，因此 (3) 式可以寫作

${\displaystyle \quad (4)\qquad \epsilon (x,t)=\sum _{m=1}^{M}e^{at}e^{ik_{m}x}}$

其中 ${\displaystyle a}$ 為常量。

由於誤差所滿足的差分方程是線性的（級數每一項的行為與整個級數一致），只估計一項的誤差變化便足以估計整體趨勢：

${\displaystyle \quad (5)\qquad \epsilon _{m}(x,t)=e^{at}e^{ik_{m}x}.}$

為找出誤差隨時間步的變化, 將方程 (5) 式應用於離散後的誤差表達式上

• ${\displaystyle {\begin{aligned}\epsilon _{j}^{n}&=e^{at}e^{ik_{m}x}\\\epsilon _{j}^{n+1}&=e^{a(t+\Delta t)}e^{ik_{m}x}\\\epsilon _{j+1}^{n}&=e^{at}e^{ik_{m}(x+\Delta x)}\\\epsilon _{j-1}^{n}&=e^{at}e^{ik_{m}(x-\Delta x)},\end{aligned}}}$

再代入到 (2) 式中，求解方程後可得

${\displaystyle \quad (6)\qquad e^{a\Delta t}=1+{\frac {\alpha \Delta t}{\Delta x^{2}}}\left(e^{ik_{m}\Delta x}+e^{-ik_{m}\Delta x}-2\right).}$

使用已知的指數三角關係式

${\displaystyle \qquad \cos(k_{m}\Delta x)={\frac {e^{ik_{m}\Delta x}+e^{-ik_{m}\Delta x}}{2}}}$ 和 ${\displaystyle \sin ^{2}{\frac {k_{m}\Delta x}{2}}={\frac {1-\cos(k_{m}\Delta x)}{2}}}$

可以將方程 (6) 變作

${\displaystyle \quad (7)\qquad e^{a\Delta t}=1-{\frac {4\alpha \Delta t}{\Delta x^{2}}}\sin ^{2}(k_{m}\Delta x/2).}$

定義漲幅因子

${\displaystyle G\equiv {\frac {\epsilon _{j}^{n+1}}{\epsilon _{j}^{n}}},}$

則誤差有限的充要條件為 ${\displaystyle \vert G\vert \leq 1}$ 。 已知

${\displaystyle \quad (8)\qquad G={\frac {e^{a(t+\Delta t)}e^{ik_{m}x}}{e^{at}e^{ik_{m}x}}}=e^{a\Delta t},}$

聯立 (7) 和 (8) 兩式，易得穩定性條件為

${\displaystyle \quad (9)\qquad \left\vert 1-{\frac {4\alpha \Delta t}{\Delta x^{2}}}\sin ^{2}(k_{m}\Delta x/2)\right\vert \leq 1}$

即

${\displaystyle \quad (10)\qquad {\frac {\alpha \Delta t}{\Delta x^{2}}}\leq {\frac {1}{2}}.}$

(10) 即為該算法的穩定性條件。 對於 FTCS 求解一維熱傳導方程，給定 ${\displaystyle \Delta x}$ ， 所允許的 ${\displaystyle \Delta t}$ 取值需要足夠小以滿足 (10) ，才能保證計算的數值穩定。