数值分析中, 冯诺依曼稳定性分析 (亦作傅立叶稳定性分析) 用于验证计算线性偏微分方程时使用特定有限差分法的数值稳定性,该分析方法基于对数值误差的傅立叶分解。1947年英国研究人员約翰·克蘭克和菲利斯·尼科爾森在文章中对该方法进行了简要介绍,
尔后又出现在冯诺依曼合作的文章中
。
洛斯阿拉莫斯国家实验室对该方法进行了进一步发展。
数值稳定性
数值稳定性与数值误差密切相关。使用有限差分方法进行计算时,若任意时间步的误差不会导致其后计算结果的发散,则可称该有限差分法是数值稳定的。如果误差随着进一步计算降低最终消失,该算法被认为稳定;若误差在进计算中保持为常量,则认为该算法“中性稳定”。但如果误差随着进一步计算增长,结果发散,则数值方法不稳定。数值方法的稳定性可以通过冯诺依曼稳定性分析得到验证。稳定性一般不易分析,特别是针对非线性偏微分方程。
冯诺依曼稳定性方法只适用于满足 Lax–Richtmyer 条件 (Lax 等价定理) 的某些特殊差分法: 偏微分方程系统须线性,常系数,满足周期性边界条件,只有两个独立变量,差分法中最多含两层时间步。 由于相对简单,人们常使用冯诺依曼稳定性分析代替其他更为详细的稳定性分析,用以估计差分方法中对容许步长的限制。
方法描述
冯诺依曼误差分析将误差分解为傅立叶级数。为了描述此过程,考虑一维热传导方程
:
\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 亦满足离散方程 :
:
\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,M, M = 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) ,才能保证计算的数值稳定。
參考資料
评论 (0)