数值分析中, 冯诺依曼稳定性分析 (亦作傅立叶稳定性分析) 用于验证计算线性偏微分方程时使用特定有限差分法的数值稳定性[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.