误差分析

第一讲误差分析¶

数值分析中的问题¶

数值计算问题可以写成如下一般化的形式：寻求解 \(x\)，满足

\[ F(x,d) = 0 \tag{1.1} \]

其中 \(d\) 是问题所依赖的参数，\(F\) 是描述解 \(x\) 和参数 \(d\) 关系的函数。

在式 (1.1) 中：

\(F\) 和 \(d\) 给定，\(x\) 待求，称为 直接问题 ；
\(F\) 和 \(x\) 给定，\(d\) 未知，称为 反问题/逆问题 ；
\(x\) 与 \(d\) 已知，\(F\) 未知，称为 识别问题 。

例：阿基米德与圆周率¶

阿基米德通过圆内接正多边形逼近圆周率，得到递推表达式：

\[ a_{n+1} = \sqrt{2 - \sqrt{4 - a_n^2}}, \quad n = 1, 2, \cdots, \quad a_1 = 1 \]

圆周率近似为 \(\pi \approx 3 \cdot 2^{n-1} a_n\)。

尽管序列 \(\{3 \cdot 2^{n-1} a_n\}\) 是 \(\pi\) 的一个很好的估计，但从数值分析的角度，仍然存在以下问题：

这个方法是在计算圆周率 \(\pi\) 吗？ (一致性问题)
序列收敛吗？ (收敛性问题)
序列 \(\{a_n\}\) 中每一项都有舍入误差，会对结果有影响吗？(例如计算机中采用 \(0.51763809020504 \ldots\) 近似 \(a_2 = \sqrt{2 - \sqrt{3}}\)) (稳定性问题)
序列要算到第 \(n\) 项时，与 \(\pi\) 误差多少？ (向前先验误差分析)
序列要算到多少项，才能保证 \(|3 \cdot 2^{n-1} a_n - \pi| < \epsilon_{tol}\)？ (向后先验误差分析)
序列要算到第 \(n\) 项时，如何根据当前的计算结果评估误差（先验方法可能无效）？ (后验误差分析)
这个方法收敛率如何？ (收敛率分析)

数值分析中的基本概念¶

问题分类¶

按问题形式分为：

直接问题：已知 \(F\) 和 \(d\)，求 \(x\)；
反问题：已知 \(F\) 和 \(x\)，求 \(d\)；
识别问题：已知 \(x\) 和 \(d\)，求 \(F\)。

适定性¶

如果定解问题的解满足以下三个条件，则称该问题是 适定的 （well-posed）：

解存在：对于给定的定解问题，至少存在一个解；
解唯一：在解存在的前提下，解是唯一的；
解连续依赖于数据：当参数发生微小变化时，解的变化也很小。

具体来说，解对参数的连续性可通过参数 \(d\) 的微扰描述：

\[ F(x + \delta x, d + \delta d) = 0 \tag{1.2} \]

满足：

\[ \forall d,\ \exists \eta_0 = \eta_0(d) > 0,\ \exists K_0 = K_0(\eta_0), \text{ s.t. } \forall \delta d: \|\delta d\| \leqslant \eta_0 \Rightarrow \|\delta x\| \leqslant K_0 \|\delta d\| \tag{1.3} \]

例题 1.1¶

求多项式 \(p(x) = x^4 - (2a - 1)x^2 + a(a - 1)\) 的实根，考察适定性。

解：该多项式在 \(a \geqslant 1\) 时有 4 个实根，在 \(1 > a \geqslant 0\) 时有 2 个实根，在 \(a < 0\) 时无实根。解的数量与参数 \(a\) 没有连续性关系，因此该问题是不适定的。

条件数¶

由式 (1.2) 可引出度量问题适定性程度的概念—— 条件数 。

定义 1.2 对于问题 (1.1) 及其扰动 (1.2)，定义：

相对条件数：

\[ \mathrm{Cond}(d) \equiv \sup_{\delta d \in \mathcal{D}} \frac{\|\delta x\| / \|x\|}{\|\delta d\| / \|d\|} \tag{1.4} \]

绝对条件数：

\[ \mathrm{Cond}_{\mathrm{abs}}(d) \equiv \sup_{\delta d \in \mathcal{D}} \frac{\|\delta x\|}{\|\delta d\|} \tag{1.5} \]

其中 \(\mathcal{D}\) 是参数 \(d\) 的容许扰动空间。

一致性¶

对于问题 (1.1)，可定义近似问题序列：

\[ F_n(x_n, d_n) = 0, \quad n \geqslant 1 \tag{1.6} \]

定义 1.3 如果有

\[ F_n(x, d) - F(x, d) \rightarrow 0 \quad \text{当 } n \rightarrow \infty \tag{1.7} \]

其中 \(x\) 是式 (1.1) 的解，则称近似问题序列 (1.6) 满足 一致性条件 。

稳定性¶

定义 1.1 对于一个数值方法，如果输入数据的误差不会造成输出数据产生较大的误差，则称该方法是稳定的；反之，若输入数据的误差会在结果中放大，则称该方法是不稳定的。

对于近似问题序列，同样可定义解的稳定性：

\[ \forall d_n,\ \exists \eta_0 = \eta_0(d_n) > 0,\ \exists K_0 = K_0(\eta_0), \text{ s.t. } \forall \delta d_n: \|\delta d_n\| \leqslant \eta_0 \Rightarrow \|\delta x_n\| \leqslant K_0 \|\delta d_n\| \tag{1.8} \]

类似地，可定义近似问题的条件数：

相对条件数：

[ \mathrm{Cond}n(d_n) \equiv \sup ]}_n} \frac{|\delta x_n| / |x_n|}{|\delta d_n| / |d_n|} \tag{1.9

绝对条件数：

[ \mathrm{Cond}{\mathrm{abs}, n}(d_n) \equiv \sup ]}_n} \frac{|\delta x_n|}{|\delta d_n|} \tag{1.10

例题 1.2¶

考察积分 \(I_n = \int_0^1 \frac{x^n}{x+5} dx\) 的数值解，分析两种方法的稳定性。

方法 1：利用递推公式 \(I_n + 5I_{n-1} = \frac{1}{n}\) 及 \(I_0 = \ln 1.2\) 逐次计算。

误差分析：\(e_n = I_n - I_n^{\ast} = -5 e_{n-1} \Rightarrow e_n = (-5)^n e_0\)
初值误差被不断放大，方法不稳定。

方法 2：利用递推公式 \(I_{k-1} = \frac{1}{5} \left[ \frac{1}{k} - I_k \right]\) 反向计算。

误差分析：\(e_{n-1} = -\frac{1}{5} e_n \Rightarrow e_0 = \left( -\frac{1}{5} \right)^n e_n\)
误差不断缩小，方法稳定。

收敛性¶

采用近似问题序列 (1.6) 求解问题 (1.1)，需研究其收敛性。

定义 1.4 对于方法 (1.6)，当且仅当满足

\[ \forall \epsilon > 0,\ \exists n_0 = n_0(\epsilon),\ \exists \delta = \delta(n_0, \epsilon), \text{ s.t. } \forall n > n_0,\ \forall \delta d_n: \|\delta d_n\| \leqslant \delta \Rightarrow \|x(d) - x_n(d + \delta d_n)\| \leqslant \epsilon \tag{1.11} \]

时，该方法是收敛的。这里 \(x(d)\) 是 (1.1) 的解，\(x_n(d + \delta d_n)\) 是 (1.6) 的解。

定理 1.1（Lax-Richtmyer 定理）

对于一个满足一致性条件的数值问题，其稳定性等价于收敛性。

先验与后验分析¶

在分析数值方法的稳定性和收敛性时，通常有以下两个策略：

向前分析（forward analysis）：通过参数 \(d\) 的扰动和数值方法的方法误差，估计解 \(x\) 的扰动大小；
向后分析（backward analysis）：通过给定 \(x\) 的目标扰动大小，基于数值方法的方法误差，估计参数 \(d\) 的容许扰动范围。

以上两个问题均属于 先验分析（priori analysis） ，通常提供在给定问题条件下数值解与真实解之间误差的上界。

后验分析（posteriori analysis） 是在数值求解之后进行的分析，侧重于利用实际计算中获得的信息来评估数值解的质量，其重要作用是进行**自适应误差控制（adaptive error control）**：

评估当前解的质量；
确定是否需要改善某些参数 \(d\) 并重新计算；
重复直至误差达到预设阈值。

误差的定义和来源¶

误差的来源¶

模型误差：实际问题的解与数学模型的解之差。常见于对物理过程进行必要的简化，例如忽略次要因素（摩擦力、空气阻力等）。
观测误差：数学问题中的参量由观测得到，但观测不可能绝对准确，由此产生的误差。
截断误差（方法误差）：由于对问题进行简化（如用有限项代替无限项）所引起的误差。例如，Taylor 展开中舍弃高阶项产生的误差。
舍入误差：计算机中的实数均为有限位小数，需通过四舍五入略去无法表示的数位，由此产生的误差。

例题 1.3¶

采用球体表面积公式 \(A = 4\pi r^2\) 计算地球表面积，取地球半径 \(r \approx 6370\ \mathrm{km}\)，分析误差来源。

解：

用球体近似地球 → 模型误差
取 \(r \approx 6370\ \mathrm{km}\) → 观测误差
\(\pi\) 取有限位、实数乘法 → 舍入误差

例题 1.4¶

对于单摆运动方程 \(L\frac{d^2\phi}{dt^2} + g\sin\phi + \mu\frac{d\phi}{dt} = 0\)，在小角度下近似为 \(L\frac{d^2\phi}{dt^2} + g\phi + \mu\frac{d\phi}{dt} = 0\)，分析误差来源。

解：

用常值 \(\mu\) 表征摩擦力 → 模型误差
测量 \(L, g, \mu\) → 观测误差
\(\sin\phi \approx \phi\) → 截断误差
数值计算 → 舍入误差

绝对误差与相对误差¶

定义 1.5 设 \(x\) 为准确值，\(x^{\ast}\) 为其近似值，称

\[ e(x^{\ast}) = x - x^{\ast} \]

为误差，也称为 绝对误差 。

注意：绝对误差 ≠ 误差的绝对值，绝对误差的符号可以为正或负。

定义 1.6 通常无法给出准确的误差值 \(e(x^{\ast})\)，但可根据测量或计算估计误差绝对值的上界 \(\epsilon\)：

\[ |e(x^{\ast})| = |x - x^{\ast}| \leqslant \epsilon \tag{1.13} \]

称 \(\epsilon\) 为近似值 \(x^{\ast}\) 的__绝对误差限__，且准确值 \(x\) 满足：

\[ x^{\ast} - \epsilon \leqslant x \leqslant x^{\ast} + \epsilon \tag{1.14} \]

或记为 \(x = x^{\ast} \pm \epsilon\)。

定义 1.7 记

\[ e_{\mathrm{T}}(x^{\ast}) = \frac{e(x^{\ast})}{x} = \frac{x - x^{\ast}}{x} \]

为近似数 \(x^{\ast}\) 的__相对误差__。

定义 1.8 由于准确值 \(x\) 未知，实际计算中也常用：

\[ e_{\mathrm{T}}(x^{\ast}) = \frac{e(x^{\ast})}{x^{\ast}} = \frac{x - x^{\ast}}{x^{\ast}} \tag{1.16} \]

作为相对误差。两者之差为 \(\frac{(x - x^{\ast})^2}{x x^{\ast}}\)，是二阶小量，可忽略。

定义 1.9 若存在正数 \(\epsilon_r\) 使得

\[ |e_{\mathrm{r}}(x^{\ast})| = |e(x^{\ast})/x^{\ast}| \leqslant \epsilon_r \tag{1.18} \]

则称 \(\epsilon_r\) 为 \(x^{\ast}\) 的__相对误差限__。注意到 \(\epsilon / |x^{\ast}|\) 是 \(x^{\ast}\) 的一个相对误差限。

基本算术过程的误差传播¶

设 \(y = f(x_1, x_2, \dots, x_n)\)，近似值为 \(x_i^{\ast}\)，则近似解为 \(y^{\ast} = f(x_1^{\ast}, x_2^{\ast}, \dots, x_n^{\ast})\)。

由微分近似，有：

绝对误差：

[ e(y^{\ast}) \approx \sum_{i=1}^n \frac{\partial f(x_1^{\ast}, \dots, x_n^{\ast})}{\partial x_i} e(x_i^{\ast}) \tag{1.22} ]

相对误差：

[ e_{\mathrm{r}}(y^{\ast}) \approx \sum_{i=1}^n \frac{\partial f}{\partial x_i} \frac{x_i^{\ast}}{y ]}} e_{\mathrm{r}}(x_i^{\ast}) \tag{1.23

基本运算的误差传播公式¶

加减法：\(f(x_1, x_2) = x_1 \pm x_2\)

[ e(x_1 \pm x_2) = e(x_1) \pm e(x_2), \quad e_{\mathrm{r}}(x_1 \pm x_2) = \frac{x_1}{x_1 \pm x_2} e_{\mathrm{r}}(x_1) \pm \frac{x_2}{x_1 \pm x_2} e_{\mathrm{r}}(x_2) ]

乘法：\(f(x_1, x_2) = x_1 x_2\)

[ e(x_1 x_2) \approx x_2 e(x_1) + x_1 e(x_2), \quad e_{\mathrm{r}}(x_1 x_2) \approx e_{\mathrm{r}}(x_1) + e_{\mathrm{r}}(x_2) ]

除法：\(f(x_1, x_2) = x_1 / x_2\)

[ e(x_1/x_2) \approx \frac{1}{x_2^2} (x_2 e(x_1) - x_1 e(x_2)), \quad e_{\mathrm{r}}(x_1/x_2) \approx e_{\mathrm{r}}(x_1) - e_{\mathrm{r}}(x_2) ]

相对误差限分析¶

加减法（\(x_1, x_2\) 同号）：

[ |e_{\mathrm{r}}(x_1 + x_2)| \leqslant \max{|e_{\mathrm{r}}(x_1)|, |e_{\mathrm{r}}(x_2)|} ]

[ |e_{\mathrm{r}}(x_1 - x_2)| \leqslant \frac{|x_1 e_{\mathrm{r}}(x_1)| + |x_2 e_{\mathrm{r}}(x_2)|}{|x_1 - x_2|} ]

乘法：

[ |e_{\mathrm{r}}(x_1 x_2)| \leqslant |e_{\mathrm{r}}(x_1)| + |e_{\mathrm{r}}(x_2)| ]

除法：

[ |e_{\mathrm{r}}(x_1/x_2)| \leqslant |e_{\mathrm{r}}(x_1)| + |e_{\mathrm{r}}(x_2)| ]

例题 1.5¶

设 \(y = x^n\)，求 \(y\) 的相对误差与 \(x\) 的相对误差之间的关系。

解：

\[ e_{\mathrm{r}}(y) \approx d(\ln x^n) = n d(\ln x) = n e_{\mathrm{r}}(x) \]

所以 \(y = x^n\) 的相对误差是 \(x\) 的相对误差的 \(n\) 倍。

有效数字¶

有效数字是近似值的一种表示方法，可同时表示近似值的大小与精确程度。

定义 1.10 如果近似值 \(x^{\ast}\) 的误差限是 \(\frac{1}{2} \times 10^{-n}\)，则称 \(x^{\ast}\) 准确到小数点后 \(n\) 位。

设 \(x^{\ast}\) 可写成规格化浮点数：

\[ x^{\ast} = \pm 0.\alpha_1\alpha_2\dots\alpha_n\dots \times 10^m \tag{1.29} \]

其中 \(\alpha_k \in \{0,1,\dots,9\}\)，\(\alpha_1 \neq 0\)。

定义 1.11 若

\[ |x - x^{\ast}| \leq \frac{1}{2} \times 10^{m - n} \tag{1.30} \]

则称 \(x^{\ast}\) 有 \(n\) 位__有效数字__，即在四舍五入意义下前 \(n\) 位都是准确的。此时绝对误差限为 \(\epsilon = \frac{1}{2} \times 10^{m - n}\)。

例题 1.6¶

求下列近似值的误差限：

\(x^{\ast} = 3587.64\) 是具有 6 位有效数字的近似值；
\(x^{\ast} = 0.0023156\) 是具有 5 位有效数字的近似值。

解：

\(x^{\ast} = 0.358764 \times 10^4\)，\(|x - x^{\ast}| \leqslant \frac{1}{2} \times 10^{4-6} = \frac{1}{2} \times 10^{-2}\)
\(x^{\ast} = 0.23156 \times 10^{-2}\)，\(|x - x^{\ast}| \leqslant \frac{1}{2} \times 10^{-2-5} = \frac{1}{2} \times 10^{-7}\)

有效数字与相对误差的关系¶

定理 1.2 若 \(x^{\ast} \neq 0\) 有 \(n\) 位有效数字，则其相对误差限满足：

\[ \frac{|x - x^{\ast}|}{|x^{\ast}|} \leqslant \epsilon_{\mathrm{r}} = \frac{1}{2\alpha_1} \times 10^{-(n-1)} \tag{1.32} \]

定理 1.3 若 \(x^{\ast} \neq 0\) 的相对误差限满足：

\[ \epsilon_{\mathrm{r}} \leqslant \frac{1}{2(\alpha_1 + 1)} \times 10^{-(n-1)} \tag{1.33} \]

则 \(x^{\ast}\) 至少有 \(n\) 位有效数字。

例题 1.7¶

用 \(x^{\ast} = 2.72\) 表示具有三位有效数字的自然对数底 \(e\) 的近似值，计算相对误差限。

解：\(\alpha_1 = 2\)，\(n = 3\)，相对误差限为 \(\frac{1}{2 \times 2} \times 10^{-(3-1)} = \frac{1}{4} \times 10^{-2}\)。

精度与准确度¶

精度：与有效数字位数有关；
准确度：与准确的有效数字位数有关。

例如，3.14111111 具有 9 位精度，但作为 \(\pi\) 的近似，只有 4 位准确的有效数字。

浮点数¶

定点数表示¶

用 \(N\) 个比特位表示实数，符号位 1 位，整数和小数部分位数固定。设小数部分占 \(k\) 位，则实数可表示为：

\[ x = (-1)^s \cdot \beta^{-k} \sum_{j=0}^{N-2} a_j \beta^j \quad (\beta = 2, s = 0,1) \tag{1.34} \]

可表示的绝对值范围：

\[ |x|_{\max} = \beta^{-k} \sum_{j=0}^{N-2} \beta^j, \quad |x|_{\min} = \beta^{-k} \tag{1.35} \]

范围非常有限。

浮点数表示¶

允许小数点位置变化，表示形式为：

\[ x = (-1)^s \cdot [0.a_1a_2\cdots a_t] \cdot \beta^e = (-1)^s \cdot m \cdot \beta^{e-t} \tag{1.36} \]

其中：

\(t\)：小数部分比特位数；
\(m\)：尾数，整数，\(0 \leqslant m \leqslant \beta^t - 1\)；
\(e\)：指数，整数，范围 \(L \leqslant e \leqslant U\)，通常 \(L < 0, U > 0\)。

为避免表示不唯一，要求 \(a_1 \neq 0\)，即 \(\beta^{t-1} \leqslant m \leqslant \beta^t - 1\)。

浮点数系统可表示的实数集合为：

\[ \mathbb{F}(\beta, t, L, U) \equiv \{0\} \bigcup \left\{ x \in \mathbb{R} \ \middle|\ x = (-1)^s \beta^e \sum_{i=1}^t a_i \beta^{-i}, a_1 \neq 0 \right\} \tag{1.37} \]

改进定义（允许 \(e = L\) 时 \(a_1 = 0\)）：

\[ \mathbb{F}_D(\beta, t, L, U) \equiv \{0\} \bigcup \{ x \in \mathbb{R} \ |\ x = (-1)^s \beta^e \sum_{i=1}^t a_i \beta^{-i}, a_1 \neq 0 \text{ if } e > L \} \]

IEEE-754 标准¶

精度	符号位 (s)	指数位 (E)	尾数位 (M)
单精度	1 位	8 位	23 位
双精度	1 位	11 位	52 位

实数表示：

单精度：\(x = (-1)^s \times (1 + M) \times 2^{E - 127}\)
双精度：\(x = (-1)^s \times (1 + M) \times 2^{E - 1023}\)

示例：1.0/3.0 的存储¶

单精度二进制：00111110101010101010101010101011
单精度实际值：0.333333343267441
双精度二进制：0011111111010101010101010101010101010101010101010101010101010101
双精度实际值：0.333333333333333

纯文本输出 vs. 二进制输出¶

纯文本输出：如 "3.33333333E-01"，占用 138 字节；
二进制输出：单精度 4 字节，双精度 8 字节，更紧凑。

在数值计算中需要注意的一些现象¶

1. 避免两个相近的数相减¶

由相对误差公式：

\[ |e_{\mathrm{r}}(x_1 - x_2)| \leqslant \frac{|x_1 e_{\mathrm{r}}(x_1)| + |x_2 e_{\mathrm{r}}(x_2)|}{|x_1 - x_2|} \quad (x_1, x_2 \text{同号}) \]

若 \(x_1\) 与 \(x_2\) 相近，分母很小，相对误差会被放大。

例题 1.8¶

当 \(x\) 充分大时，计算 \(\sqrt{x+1} - \sqrt{x}\) 应改为： [ \sqrt{x+1} - \sqrt{x} = \frac{1}{\sqrt{x+1} + \sqrt{x}} ] 以减少相对误差。

2. 防止数量级相差很大的数进行加减运算¶

由于浮点数尾数位数有限，数量级相差太大的数相加时，绝对值较小的数可能失效。

例题 1.9¶

假设计算机只能存储 8 位有效数字，用求根公式计算方程 \(x^2 - (10^9 + 1)x + 10^9 = 0\) 的根（理论解 \(x_1 = 10^9, x_2 = 1\)）。

解：计算 \(b^2 - 4ac\) 时，\(b^2\) 与 \(4ac\) 相差极大，导致 \(b^2 - 4ac \approx b^2\)，从而 \(x_2 \approx 0\)。

改进算法：

\[ x_1 = \frac{-b - \mathrm{sign}(b)\sqrt{b^2 - 4ac}}{2a}, \quad x_2 = \frac{c}{a x_1} \]

3. 注意简化计算步骤¶

简化步骤可降低计算复杂度，避免误差累积。

例题 1.10¶

计算多项式 \(P(x) = a_n x^n + a_{n-1} x^{n-1} + \dots + a_0\)。

解：依次计算需 \(\frac{n(n+1)}{2}\) 次乘法与 \(n\) 次加法。 秦九韶算法 ：

\[ P(x) = (\dots((a_n x + a_{n-1})x + a_{n-2})x + \dots)x + a_0 \]

只需 \(n\) 次乘法与 \(n\) 次加法。

4. 选择合适的精度类型¶

阿丽亚娜5型火箭失事（1996年）：因 64 位浮点数转换为 16 位有符号整数时溢出，导致惯性导航系统故障，火箭坠毁。
爱国者导弹拦截失败（1991年）：计算机开机后时间计算的累积误差导致拦截失败。

5. 重视累积误差¶

长时间运行的计算中，微小误差可能累积到不可接受的程度，需通过算法设计或定期校正来控制。

第一讲 误差分析¶