1.2 误差知识与算法知识

本节是全书的地基，只讲一件事：误差是怎么走路的。具体分四步：误差从哪里来（1.2.1）；怎么度量一个近似值好不好（1.2.2）；它如何跟着运算一步步传下去、被放大（1.2.3–1.2.4）；以及怎么设计算法把它管住（1.2.4–1.2.6）。原书例 1–例 5 与全部经典论证完整保留，证明与推导全部展开。

1.2.1 误差的来源与分类

工程计算中，估计结果的精确度十分重要，而精确度受各种误差的影响。按来源划分，误差有四类。

典型例子：Taylor 展开的截断。由 Taylor 公式，函数 $f(x)$ 在 $x=0$ 处可表示为

$$f(x)=f(0)+f'(0)x+\frac{f''(0)}{2!}x^2+\cdots+\frac{f^{(n)}(0)}{n!}x^n+\frac{f^{(n+1)}(\theta x)}{(n+1)!}x^{n+1},\qquad 0\lt\theta\lt1 .$$

当 $\lvert x\rvert$ 不大时，去掉右端最后一项（余项），得近似公式

$$f(x)\approx f(0)+f'(0)x+\frac{f''(0)}{2!}x^2+\cdots+\frac{f^{(n)}(0)}{n!}x^n .$$

此近似公式的误差就是截断误差——它来自「方法」对无穷级数尾部的舍弃，与机器无关。例如 $\sin x \approx x - x^3/6$ 在 $x=0.5$ 处的截断误差约为 $0.5^5/120 \approx 2.6\times10^{-4}$。

原书给出两个十位十进制下的例子。第一个是 $1\div3 = 0.3333333333$：无穷循环小数被截断。

第二个是 $(1.000002)^2-1.000004 = 0$，而精确值是 $4\times10^{-12}$：$(1.000002)^2 = 1.000004000004$，减去 $1.000004$ 应剩 $0.000000000004$；但在十位字长下，两个操作数都只保留到第 10 位数字，相减结果被舍成 $0$。这里产生的误差就是舍入误差。

1.2.2 绝对误差、相对误差与有效数字

绝对误差与误差限

准确值 $x$、近似值 $a$ 和误差限 $\varepsilon$ 三者的关系是

$$a-\varepsilon\lt x\lt a+\varepsilon ,$$

常记为 $x = a\pm\varepsilon$。例如，$a=3.14$ 作为圆周率 $\pi$ 的一个近似值，绝对误差为 $e=\pi-3.14$。由 $\pi=3.14159265\cdots$ 易知 $\lvert e\rvert\lt0.002$，故可取误差限 $\varepsilon = 0.002$，写作 $\pi = 3.14\pm0.002$。

用绝对误差刻画精确程度有局限：它没有反映误差在数本身中所占的比例。用最小刻度为毫米的卡尺量一根火柴（$\approx20$ mm）和量一根竹竿（$\approx2$ m），绝对误差限同样是 $0.5$ mm，但两者的可信程度天差地别——这就需要相对误差。

相对误差与相对误差限

两种相对误差差多远？（推导）真实值 $(x-a)/x$ 与实用值 $(x-a)/a$ 之间只差一个因子：

$$\frac{x-a}{a}=\frac{x-a}{x}\cdot\frac{x}{a}=e_{\mathrm r}\cdot\frac{1}{1-e_{\mathrm r}} = e_{\mathrm r},(1+e_{\mathrm r}+e_{\mathrm r}^2+\cdots),$$

可见二者之差是 $e_{\mathrm r}$ 的高阶小量 $e_{\mathrm r}^2+e_{\mathrm r}^3+\cdots$。数值计算里 $\lvert e_{\mathrm r}\rvert$ 通常远小于 $1$（否则近似值毫无意义），替换是安全的；只有 $\lvert e_{\mathrm r}\rvert\ge0.1$ 量级时才需要警惕，此时近似值本身已接近失效。另外，$x=0$ 处相对误差无定义——近似 $0$ 只能靠绝对误差谈论。

有效数字

三条直接推论（原书所述，逐条解释）：

1. 由准确值经四舍五入得到的近似值，从它的末位数字到第一位非零数字都是有效数字（例 2 的结论换一种说法）。
2. 同一个准确值的不同近似值，有效数字越多，绝对误差与相对误差都越小。绝对误差限 $\tfrac12\times10^{m-n}$ 随 $n$ 增大而按十倍速缩小；相对误差限见下面的 定理 1.2.1，也随 $n$ 缩小。一句话版本：每多一位有效数字，相对误差缩小十倍（3 位 ⟹ $\lesssim0.5$，6 位 ⟹ $\lesssim0.0005$）。
3. 尾零有意义：$2.140012$ 的两个近似值 $2.14$ 与 $2.1400$ 写法有区别——前者 3 位有效数字，后者 5 位。准确值的有效数字可视为有无限多位。

1.2.3 函数求值的误差估计

设 $u=f(x)$ 存在足够高阶的导数，$a$ 是自变量 $x$ 的近似值，则 $\tilde u=f(a)$ 是函数值 $u=f(x)$ 的近似值。本小节回答：自变量的误差 $e(a)$ 如何决定函数值的误差 $e(\tilde u)$？

一元情形

放大因子怎么读。$\lvert f'(a)\rvert$ 衡量误差从输入传到输出被放大多少倍：$\lvert f'(a)\rvert\le1$ 时，函数反而「吸收」误差；$\lvert f'(a)\rvert\gg1$ 时，输入的小误差被放大成输出的大误差。这正是「问题的条件数」概念的雏形，1.2.6 会把它精确化。

多元情形

设 $n$ 元函数 $u=f(x_1,x_2,\cdots,x_n)$ 充分可微，$a_i$ 是 $x_i$ 的近似值，则 $\tilde u=f(a_1,a_2,\cdots,a_n)$ 是 $u$ 的近似值。由多元函数 Taylor 公式（一阶展开）得

$$e(\tilde u)\approx\sum_{i=1}^{n}\frac{\partial f(a_1,a_2,\cdots,a_n)}{\partial x_i},e(a_i),\qquad \varepsilon(\tilde u)\approx\sum_{i=1}^{n}\left\lvert\frac{\partial f(a_1,a_2,\cdots,a_n)}{\partial x_i}\right\rvert\varepsilon(a_i). \tag{1.1}$$

推导要点：把其余变量视为常数、对每个变量分别用一元 Taylor 展开（一阶），相加后二阶交叉项 $O([e(a_i)]^2)$ 与 $O(e(a_i)e(a_j))$ 在 $\lvert e(a_i)\rvert$ 都很小时可略。

如果所有 $\left\lvert\partial f/\partial x_i\right\rvert$ 全为零或全都很小（一阶项消失），则要使用 Taylor 公式中的高阶项——典型场景正是下面四则运算中的「相近数相减」。

四则运算的误差估计

设 $a$ 和 $b$ 分别是准确值 $x$ 和 $y$ 的近似值，则 $a+b,\ a-b,\ ab,\ a/b\ (b\neq0)$ 分别是 $x+y,\ x-y,\ xy,\ x/y$ 的近似值。对四则运算逐个使用式 (1.1)：

1.2.4 算法及其计算复杂性

先把「算法」这个词说准，再谈用什么标准衡量算法的优劣。

衡量算法优劣有两个标准：一是要有可靠的理论基础，包括正确性、收敛性、数值稳定性，并可作误差分析；二是要有良好的计算复杂性。所谓计算复杂性，是指达到给定精度时该算法所需的计算量（时间复杂性）与所占的内存空间（空间复杂性）。

例行公事：多项式求值的两种算法

计算 $p(x)=\displaystyle\sum_{j=0}^{n}a_jx^j$ 的值，输入 $a_j\ (j=0,1,\cdots,n)$ 和 $x$，输出 $p(x)$。

实测对照（$n=5000$，随机系数，$x=0.9999$，本教材真实运行数据）：

import time, numpy as np
rng = np.random.default_rng(0)
a = rng.uniform(-1, 1, 5001)          # a[k] 是 x^k 的系数
x = 0.9999

t0 = time.perf_counter()              # 算法一：幂每次从 1 重乘
s = 0.0
for k in range(5001):
    xk = 1.0
    for _ in range(k):
        xk = xk * x
    s += a[k] * xk
t1 = time.perf_counter()

T = 0.0                               # 算法二：秦九韶
for k in range(5000, -1, -1):
    T = T * x + a[k]
t2 = time.perf_counter()

print(f"算法一: {(t1-t0)*1000:.0f} ms   算法二: {(t2-t1)*1000:.2f} ms")

算法一: 279 ms   算法二: 0.39 ms

约 700 倍。理论乘法次数之比 $n(n+1)/2 : n = (n+1)/2 : 1 \approx 2500:1$，实测的 700 倍与它同数量级（常数因子来自 Python 解释器开销）。规模越大差距越大：$O(n^2)$ 与 $O(n)$ 的差别是结构性的。

设计算法的四条原则

算法通常在计算机上执行，而计算机的字长有限，舍入误差因而不可避免。为了减少舍入误差的影响，设计算法时应遵循以下原则。

原则 (1) 要有数值稳定性，即能控制舍入误差的传播

原书经典反例：一个递推的正反两个方向。要在四位十进制的限制下计算积分

$$y_n=\int_0^1\frac{x^n}{x+5},dx\qquad(n=0,1,\cdots,100).$$

建立递推关系（推导）。由 $x^n+5x^{n-1}=x^{n-1}(x+5)$，分子恰好被分母整除：

$$y_n+5y_{n-1}=\int_0^1\frac{x^n+5x^{n-1}}{x+5},dx=\int_0^1 x^{n-1},dx=\frac1n\quad(n\ge1),$$

$$y_0=\int_0^1\frac{dx}{x+5}=\ln 6-\ln 5\approx0.1823 .$$

正向递推（不稳定）：$y_n=\dfrac1n-5y_{n-1}$。误差分析：设 $y_{n-1}$ 带误差 $\delta_{n-1}$，则算出的 $y_n$ 带误差 $\delta_n=-5\delta_{n-1}$——每步放大 5 倍。初值 $y_0=0.1823$ 的舍入误差约 $-2.2\times10^{-5}$（真值 $0.18232156\cdots$），传到 $y_{100}$ 时放大 $5^{100}\approx7.9\times10^{69}$ 倍。数值实验（双精度执行，递推系数是元凶，与字长无关）：

 n        前推值                真值             误差
 0   1.823000e-01      1.823216e-01     -2.2e-05     ← 初值只差 5 位
 4   2.083333e-02      3.430633e-02     -1.3e-02
 6  -3.125000e-01      2.432491e-02     -3.4e-01     ← 已失真
10  -2.105002e+02      1.536755e-02     -2.1e+02
30  -2.007633e+16      5.404633e-03     -2.0e+16     ← 荒谬

先给反向递推找一个起点。反向递推要从 $y_{100}$ 出发，可它的真值并不知道——好在真值的范围可以估出来。在 $[0,1]$ 上 $5\le x+5\le6$，故 $x^n/6\le x^n/(x+5)\le x^n/5$，积分得

$$\frac{1}{6(n+1)}\lt y_n\lt\frac{1}{5(n+1)} .$$

取 $y_{100}\approx\dfrac12\left(\dfrac{1}{606}+\dfrac{1}{505}\right)\approx0.001815$（两端中点，误差不超过区间宽度之半 $\approx1.7\times10^{-4}$）。

反向递推（稳定）：把递推关系反过来，$y_{n-1}=\dfrac{1}{5n}-\dfrac{1}{5}y_n$：

$$\begin{cases}y_{100}\approx0.001815,\ y_{n-1}=\dfrac{1}{5n}-\dfrac15 y_n\quad(n=100,\ 99,\ \cdots,\ 1).\end{cases}$$

每步误差缩小为 1/5：$\lvert\delta_{n-1}\rvert=\lvert\delta_n\rvert/5$。起点误差 $1.7\times10^{-4}$ 传到 $y_0$ 时只剩 $1.7\times10^{-4}/5^{100}$，加上途中新产生的舍入误差（同样每步被 5 压制），实测 $y_0$ 的误差为 $-2.8\times10^{-17}$——机器精度级别。

原则 (2) 两数相加时，防止较小的数加不到较大的数中

较小数加不到较大数中有时是允许的，但有时会产生严重后果。原书例：在十位十进制的限制下求解一元二次方程

$$x^2+10^4x-0.01=0 ,$$

使用求根公式 $x=\dfrac{-b\pm\sqrt{b^2-4ac}}{2a}$。按加法的对阶（alignment）规则：

$$b^2-4ac=10^8+0.04=0.1\times10^9+0.00000000004\times10^9 ,$$

十位字长下 $0.00000000004\times10^9$ 被当做 $0$，故 $b^2-4ac=10^8$，于是解得

$$x_1=0,\qquad x_2=-10^4 .$$

$x_2=-10^4$ 是合理的、可接受的；但 $x_1=0$ 不可接受（真值 $10^{-6}$，误差 100%），后果还很严重：下游计算只要拿 $x_1$ 作除数就会崩溃。补救的办法是利用根与系数的关系 $x_1x_2=c/a=-0.01$：

$$x_1=-\frac{0.01}{x_2}=\frac{0.01}{10^4}=10^{-6} .$$

现代版本：同样的灾难在 IEEE 双精度里照样发生。方程 $x^2-10^8x+1=0$（两根 $10^{\pm8}$）用 float64 直接套求根公式，小根算得 $x_2=7.45\times10^{-9}$，相对误差 $25.5$——病根是判别式开方时 $10^8-\sqrt{10^{16}-4}$ 发生相近数相减；float32 下小根干脆得 $0.0$。改用韦达公式 $x_2=c/(a,x_1)$ 重算小根，float64 与 float32 都精确得到 $10^{-8}$（相对误差 $0$）。完整实验见 1.2.7 实验二。

原则 (3) 尽量避免两个相近的近似值相减，以免严重损失有效数字

原书例：$x=1.232,\ y=1.231$ 是两个准确值，在四位十进制限制下计算 $z=x^3-y^3$。

直接算法：$z=1.232^3-1.231^3\approx1.870-1.865=0.005$。真值 $z=0.004550$，误差 $4.5\times10^{-4}$ 达到结果本身的 $10$——最多一位有效数字。原因：$1.870$ 与 $1.865$ 是相近数相减（各自带 $\pm0.0005$ 的舍入，相减后 $0.005$ 中混入了最多 $0.001$ 的误差）。

稳定算法：因式分解 $x^3-y^3=(x-y)(x^2+xy+y^2)$，其中 $x-y=0.001$ 是准确值：

$$z\approx0.001\times(1.518+1.517+1.515)=0.004550 .$$

新误差只来自三个平方 / 乘积项的舍入（各 $\pm0.0005$）：

$$\varepsilon(z)=0.001\times(0.0005+0.0005+0.0005)=0.0000015\lt0.000005 ,$$

故 $z=0.004550$ 至少有三位有效数字——从 1 位到 3 位，代价只是一次代数变形。

一个重要的现代澄清（Sterbenz 引理）：