数值稳定性

在数值分析中，数值稳定性是一种希望得到的数值算法特性。根据算法的不同，稳定性的精确定义也有所不同，但是都与算法的精确性与正确性相关。

理论上有些计算下可以用多种代数上等价的理想实数或者复数算法来实现，但是实际上由于不同的数值稳定性可能会得到不同的结果。数值稳定性的一项任务就是选择強健（robust，即有良好数值稳定性）的算法。

在数值线性代数中经常使用前向、后向以及混合稳定性的概念。

假设要用数值算法解决的问题是用函数${\displaystyle f}$将数据${\displaystyle x}$映射到解${\displaystyle y}$，通常算法的结果${\displaystyle y^{*}}$会与“真”解${\displaystyle y}$有一定的偏差。误差的来源主要有舍入误差、截断误差以及数据误差。算法的「前向误差」是结果与真解之间的差别，即${\displaystyle \Delta y=y^{*}-y}$。「后向误差」是满足${\displaystyle f(x+\Delta x)=y^{*}}$的最小${\displaystyle \Delta x}$，也就是说后向误差说明算法的所解决的问题。前向误差和后向误差通过条件数发生关系：前向误差的幅度最多是条件数乘以后向误差的幅度。

在许多情况下，需要考虑相对误差

${\displaystyle {\frac {|\Delta x|}{|x|}},}$

而不是绝对误差${\displaystyle \Delta x}$。

如果对于任意的输入${\displaystyle x}$来说后向误差都很小，那么算法就是「后向稳定」的。当然，“小”是一个相对的概念，需要根据所用的场合进行定义。通常要求误差要与单位舍入误差（英语：Machine epsilon）处于同一数量级，或是略大一些。

通常数值稳定性的定义是使用一个包括了前向误差与后向误差的更加宽泛的概念，称为「混合稳定性」。按照这个概念，如果一个算法是稳定的，那么存在${\displaystyle \Delta x}$使得${\displaystyle \Delta x}$与${\displaystyle f(x+\Delta x)-y^{*}}$都很小。因此，后向稳定算法永远是稳定的。

如果算法的前向误差除以条件数得到的结果很小，那么这个算法就是「前向稳定」的。这就意味着如果一个算法的前向误差与后向稳定算法的误差幅度类似那么就是前向稳定的。

上面的定义在截断误差不重要的情况下是很确切的。在微分方程等另外一些场合中，则需要另外的数值稳定性定义。

在数值常微分方程中，有不同的数值稳定性概念，如A稳定性等。它们通常与动力系统中的李雅普诺夫稳定性等稳定性概念相关。在解刚性方程的时候稳定方法的使用很重要。

在数值偏微分方程中有另外一种数值稳定性的定义。如果随着步长逐渐趋近于零，偏微分方程的数值解仍然保持有界，那么这个算法就是稳定的。拉克斯等价定理表明如果算法是一致的，且是稳定的，那么这个算法就会收敛。有时候将数值扩散考虑在内来实现稳定性。数值扩散是一个数学术语，它保证舍入误差以及其它误差在计算的过程中逐渐散去，而不会累积起来越变越大。

根號2的計算（大約是1.41421）是適定性問題。許多演算法會用起始的近似值x₀開始，最後收斂到${\displaystyle {\sqrt {2}}}$，例如選x₀ = 1.4，接著薒計算x₁, x₂等。其中一個方法是著名的巴比倫法（英语：Babylonian method），公式為x_k+1 = (x_k+ 2/x_k)/2。另一個方法，先稱為X方法，公式是x_k+1 = (x_k² − 2)² + x_k。^{[註 1]}以下是一些迭代後的結果，初始值x₀分別是1.4和1.42。

不論初始值為何，巴比倫法都可以快速收斂。而X方法在初始值x₀ = 1.4時收斂，但收斂很慢，而在初始值為1.42時會發散。因此巴比倫法具有數值穩定的特性，而X方法沒有。

數值穩定性也會受計算時有效位元的個數所影響。假設有電腦可以保持四位的有效數字，以下二個定效的函數可以說明有效位丟失（loss of significance）的問題

${\displaystyle f(x)=x\left({\sqrt {x+1}}-{\sqrt {x}}\right)}$和${\displaystyle g(x)={\frac {x}{{\sqrt {x+1}}+{\sqrt {x}}}}.}$

比較以下兩個結果

${\displaystyle f(500)=500\left({\sqrt {501}}-{\sqrt {500}}\right)=500\left(22.38-22.36\right)=500(0.02)=10}$

和

${\displaystyle {\begin{alignedat}{3}g(500)&={\frac {500}{{\sqrt {501}}+{\sqrt {500}}}}\\&={\frac {500}{22.38+22.36}}\\&={\frac {500}{44.74}}=11.17\end{alignedat}}}$

比較上述二個結果，可以看到第一式的有效位丟失（原因是來自兩個相近數字${\displaystyle {\sqrt {501}}}$和${\displaystyle {\sqrt {500}}}$近似值相減後帶來的灾难性抵消），這對結果造成很大的影響，而這兩個函數其實是等效的

${\displaystyle {\begin{alignedat}{4}f(x)&=x\left({\sqrt {x+1}}-{\sqrt {x}}\right)\\&=x\left({\sqrt {x+1}}-{\sqrt {x}}\right){\frac {{\sqrt {x+1}}+{\sqrt {x}}}{{\sqrt {x+1}}+{\sqrt {x}}}}\\&=x{\frac {({\sqrt {x+1}})^{2}-({\sqrt {x}})^{2}}{{\sqrt {x+1}}+{\sqrt {x}}}}\\&=x{\frac {x+1-x}{{\sqrt {x+1}}+{\sqrt {x}}}}\\&=x{\frac {1}{{\sqrt {x+1}}+{\sqrt {x}}}}\\&={\frac {x}{{\sqrt {x+1}}+{\sqrt {x}}}}\\&=g(x)\end{alignedat}}}$

若用無限精度計算，其數值是11.174755...^{[註 2]}。

• 計算變異的演算法（英语：Algorithms for calculating variance）
• 穩定性理論
• 混沌理论
• 误差传播

1. ^ 這是針對${\displaystyle x=(x^{2}-2)^{2}+x=f(x)}$方程的定點迭代（英语：fixed point iteration），其解包括${\displaystyle {\sqrt {2}}}$。因${\displaystyle f(x)\geq x}$，此迭代會一直增加，在${\displaystyle x_{1}=1.4<{\sqrt {2}}}$時，會收斂，在${\displaystyle x_{1}=1.42>{\sqrt {2}}}$時則會發散。
2. ^ 這是源自Mathews & Fink (1999)修改而來^[1]

• Nicholas J. Higham, Accuracy and Stability of Numerical Algorithms, Society of Industrial and Applied Mathematics, Philadelphia, 1996. ISBN 0-89871-355-2.