冯诺依曼稳定性分析

维基百科,自由的百科全书
跳转至: 导航搜索

数值分析中, 冯诺依曼稳定性分析 (亦作傅立叶稳定性分析) 用于验证计算线性偏微分方程时使用特定有限差分法数值稳定性[1],该分析方法基于对数值误差的傅立叶分解。1947年英国研究人员 John CrankPhyllis Nicolson 在文章中对该方法进行了简要介绍[2], 尔后又出现在冯诺依曼合作的文章中 [3]洛斯阿拉莫斯国家实验室对该方法进行了进一步发展。

数值稳定性[编辑]

数值稳定性与数值误差密切相关。使用有限差分方法进行计算时,若任意时间步的误差不会导致其后计算结果的发散,则可称该有限差分法是数值稳定的。如果误差随着进一步计算降低最终消失,该算法被认为稳定;若误差在进计算中保持为常量,则认为该算法“中性稳定”。但如果误差随着进一步计算增长,结果发散,则数值方法不稳定。数值方法的稳定性可以通过冯诺依曼稳定性分析得到验证。稳定性一般不易分析,特别是针对非线性偏微分方程。

冯诺依曼稳定性方法只适用于满足 Lax–Richtmyer 条件 (Lax 等价定理) 的某些特殊差分法: 偏微分方程系统须线性,常系数,满足周期性边界条件,只有两个独立变量,差分法中最多含两层时间步[4]。 由于相对简单,人们常使用冯诺依曼稳定性分析代替其他更为详细的稳定性分析,用以估计差分方法中对容许步长的限制。

方法描述[编辑]

冯诺依曼误差分析将误差分解为傅立叶级数。为了描述此过程,考虑一维热传导方程


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

空间网格间隔为 L, 对网格作 FTCS (Forward-Time Central-Space,时间步前向欧拉法,空间步三点中心差分) 离散处理,


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

其中 r = \frac{\alpha\, \Delta t}{\Delta x^2}u_j^{n} 为离散网格上的数值解,用于近似此偏微分方程的精确解 u(x,t)

定义舍入误差 \epsilon_j^n = N_j^n - u_j^n 。 其中 u_j^n 是离散方程 (1) 式的精确解,N_j^n 为包含有限浮点精度的数值解。 因为精确解 u_j^n 满足离散方程, 误差 \epsilon_j^n 亦满足离散方程 [5]


  \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) 中,误差和数值解随时间具有一致的变化趋势。对于含周期性边界条件的线性微分方程,间隔 L 上的空间部分误差可展开为傅立叶级数


  \quad (3) \qquad \epsilon(x) = \sum_{m=1}^{M} A_m e^{ik_m x}

其中波数 k_m = \frac{\pi m}{L}m = 1,2,\ldots,MM = L/\Delta x。 通过假设误差幅度 A_m 是时间的函数,可以给出误差和时间的关系。 不难知单步中,误差随时间指数增长,因此 (3) 式可以写作


  \quad (4) \qquad \epsilon(x,t) = \sum_{m=1}^{M} e^{at} e^{ik_m x}

其中 a 为常量。

由于误差所满足的差分方程是线性的(级数每一项的行为与整个级数一致),只估计一项的误差变化便足以估计整体趋势:


  \quad (5) \qquad \epsilon_m(x,t) = e^{at} e^{ik_m x}.

为找出误差随时间步的变化, 将方程 (5) 式应用于离散后的误差表达式上

  • 
\begin{align}
 \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{align}

再代入到 (2) 式中,求解方程后可得


  \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).

使用已知的指数三角关系式

 \qquad \cos(k_m \Delta x) = \frac{e^{ik_m \Delta x} + e^{-ik_m \Delta x}}{2}  \sin^2\frac{k_m \Delta x}{2} = \frac{1 - \cos(k_m \Delta x)}{2}

可以将方程 (6) 变作


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

定义涨幅因子


  G \equiv \frac{\epsilon_j^{n+1}}{\epsilon_j^n},

则误差有限的充要条件为 \vert G \vert \leq 1 。 已知


  \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) 两式,易得稳定性条件为


  \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


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

(10) 即为该算法的稳定性条件。 对于 FTCS 求解一维热传导方程,给定 \Delta x , 所允许的 \Delta t 取值需要足够小以满足 (10) ,才能保证计算的数值稳定。

參考資料[编辑]

  1. ^ Analysis of Numerical Methods by E. Isaacson, H. B. Keller
  2. ^ 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 
  3. ^ Charney, J. G.; Fjørtoft, R.; von Neumann, J., Numerical Integration of the Barotropic Vorticity Equation, Tellus. 1950, 2: 237–254 
  4. ^ Smith, G. D., Numerical Solution of Partial Differential Equations: Finite Difference Methods, 3rd ed.. 1985:  67–68 
  5. ^ Anderson, J. D., Jr.. Computational Fluid Dynamics: The Basics with Applications. McGraw Hill. 1994.