數值分析中, 馮諾依曼穩定性分析 (亦作傅立葉穩定性分析) 用於驗證計算線性偏微分方程時使用特定有限差分法的數值穩定性[1],該分析方法基於對數值誤差的傅立葉分解。1947年英國研究人員約翰·克蘭克和菲利斯·尼科爾森在文章中對該方法進行了簡要介紹[2],
爾後又出現在馮諾依曼合作的文章中
[3] 。
洛斯阿拉莫斯國家實驗室對該方法進行了進一步發展。
數值穩定性
數值穩定性與數值誤差密切相關。使用有限差分方法進行計算時,若任意時間步的誤差不會導致其後計算結果的發散,則可稱該有限差分法是數值穩定的。如果誤差隨着進一步計算降低最終消失,該算法被認為穩定;若誤差在進計算中保持為常量,則認為該算法「中性穩定」。但如果誤差隨着進一步計算增長,結果發散,則數值方法不穩定。數值方法的穩定性可以通過馮諾依曼穩定性分析得到驗證。穩定性一般不易分析,特別是針對非線性偏微分方程。
馮諾依曼穩定性方法只適用於滿足 Lax–Richtmyer 條件 (Lax 等價定理) 的某些特殊差分法: 偏微分方程系統須線性,常係數,滿足周期性邊界條件,只有兩個獨立變量,差分法中最多含兩層時間步[4]。 由於相對簡單,人們常使用馮諾依曼穩定性分析代替其他更為詳細的穩定性分析,用以估計差分方法中對容許步長的限制。
方法描述
馮諾依曼誤差分析將誤差分解為傅立葉級數。為了描述此過程,考慮一維熱傳導方程
![{\displaystyle {\frac {\partial u}{\partial t}}=\alpha {\frac {\partial ^{2}u}{\partial x^{2}}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1aec8a435e5ce3588f59f31116e9c92687faee4b)
空間網格間隔為
, 對網格作 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)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/11cc9be387bc139fe3812888dc270134dfe784e2)
其中
。
為離散網格上的數值解,用於近似此偏微分方程的精確解
。
定義捨入誤差
。
其中
是離散方程 (1) 式的精確解,
為包含有限浮點精度的數值解。 因為精確解
滿足離散方程, 誤差
亦滿足離散方程 [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)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/5d15fc1b26535463cf09d8914da470cf60c84d39)
此式將確定誤差的遞推關係。方程 (1) 和 (2) 中,誤差和數值解隨時間具有一致的變化趨勢。對於含周期性邊界條件的線性微分方程,間隔
上的空間部分誤差可展開為傅立葉級數
![{\displaystyle \quad (3)\qquad \epsilon (x)=\sum _{m=1}^{M}A_{m}e^{ik_{m}x}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c0aac6f8bc5937ca96ae6e6dd89b589eb302eb2f)
其中波數
,
,
。 通過假設誤差幅度
是時間的函數,可以給出誤差和時間的關係。 不難知單步中,誤差隨時間指數增長,因此 (3) 式可以寫作
![{\displaystyle \quad (4)\qquad \epsilon (x,t)=\sum _{m=1}^{M}e^{at}e^{ik_{m}x}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d506c67ab59d4df1099c39836b528dd9162eeac5)
其中
為常量。
由於誤差所滿足的差分方程是線性的(級數每一項的行為與整個級數一致),只估計一項的誤差變化便足以估計整體趨勢:
![{\displaystyle \quad (5)\qquad \epsilon _{m}(x,t)=e^{at}e^{ik_{m}x}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e677b69a8659a4b1fe9a9b40e846ff7f5e3bab43)
為找出誤差隨時間步的變化, 將方程 (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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b37aa411765cdd678a0fb88f222b43c24e969508)
再代入到 (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).}](https://wikimedia.org/api/rest_v1/media/math/render/svg/eacd5ec18e35d60cd0c67ba1b392f8888c3b99da)
使用已知的指數三角關係式
和 ![{\displaystyle \sin ^{2}{\frac {k_{m}\Delta x}{2}}={\frac {1-\cos(k_{m}\Delta x)}{2}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/bb9321ea92d733051668fc072e4c2aa349674c4b)
可以將方程 (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).}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2f5aa9db791629d6388661d020b37651dbf7d665)
定義漲幅因子
![{\displaystyle G\equiv {\frac {\epsilon _{j}^{n+1}}{\epsilon _{j}^{n}}},}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a5969045388c66f1fc5d24acfb7957b25ccea32c)
則誤差有限的充要條件為
。 已知
![{\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},}](https://wikimedia.org/api/rest_v1/media/math/render/svg/915218ff007443875a8ec68f29ce8ec00976edd0)
聯立 (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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/853f72b12200f278056ae94f723fb0583aba71ce)
即
![{\displaystyle \quad (10)\qquad {\frac {\alpha \Delta t}{\Delta x^{2}}}\leq {\frac {1}{2}}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/cf1e2d036c5fb6ece08a81fd8c6f216b0d01df1f)
(10) 即為該算法的穩定性條件。 對於 FTCS 求解一維熱傳導方程,給定
, 所允許的
取值需要足夠小以滿足 (10) ,才能保證計算的數值穩定。
參考資料
- ^ Analysis of Numerical Methods by E. Isaacson, H. B. Keller. [2011-05-20]. (原始內容存檔於2011-05-21).
- ^
Crank, J.; Nicolson, P., A Practical Method for Numerical Evaluation of Solutions of Partial Differential Equations of Heat Conduction Type, Proc. Camb. Phil. Soc., 1947, 43: 50–67, doi:10.1007/BF02127704
- ^
Charney, J. G.; Fjørtoft, R.; von Neumann, J., Numerical Integration of the Barotropic Vorticity Equation, Tellus, 1950, 2: 237–254
- ^
Smith, G. D., Numerical Solution of Partial Differential Equations: Finite Difference Methods, 3rd ed.: 67–68, 1985
- ^ Anderson, J. D., Jr. Computational Fluid Dynamics: The Basics with Applications. McGraw Hill. 1994.