Lecture 1: Introduction¶
误差¶
例子：估计 \(\int_0^1 e^{-x^2} \mathrm x\)¶
我们可以使用泰勒展开，来计算这个公式：
在这里，有两个误差，一个是拉格朗日余项，一个是数值精度。
对于拉格朗日余项（\(f(x) = e^{-x}\)）：
从而：
对于数值精度，由于我们采用 rounding，因此每个存在精度损失的数值的精度损失不超过 \(0.0005\)。由于有精度损失的数值有两个：\(1/3\)，\(1/42\)，因此最大损失为 \(0.001\)。
从而，总体最大误差，至多为 \(0.006\)。
注意¶
为何选择三位精度？因为 Lagrange 余项的损失上限为 0.005。我们选择的精度的精度损失上限应该和 Lagrange 余项损失上限相符。
误差来源¶
误差来源，有三种
- rounding error
- 3.97 -> 3, 4.01 -> 4
- chopping error
- 3.49 -> 3, 3.51 -> 4
- 数据本身的误差
绝对误差和相对误差¶
相对误差在同级相减时容易放大，绝对误差在大除小时容易放大。
相对误差¶
例子：有效位数¶
对于 \(k\) 位有效位数的数字
- 如果采用 chopping，\(p\) 的相对精度是：\(10^{-k+1}\)
- e.g. \(k=2: 10.99 \to 10\)，相对精度就是 \(0.099 \approx 0.1 = 10^{-2 + 1}\)
- 如果采用 rounding，\(p\) 的相对精度是：\(0.5 * 10^{-k+1}\)
- e.g. \(k=2: 10.49 \to 10\)，相对精度就是 \(0.049 \approx 0.05 = 0.5 * 10^{-2 + 1}\)
例子：同级相减¶
对于差不多大的数字，相减可能会导致有效数字减少。以下是一个极端的例子：
- \(a = 3.1415925, \bar a = 3.1415925456789\dots\)
- \(b = 3.1415926, \bar b = 3.141592556789\dots\)
- 那么，\(\lvert a - b \rvert = 0.0000001, \lvert \bar a - \bar b \rvert = -0.000000011110\dots\)
- 从而，相对误差从之前的 \(10^{-7}\) 量级，增大到了 \(10^{-1}\) 之巨
这就是因为：
$$ \begin{aligned} &a = \bar a + e_a, b = \bar b + e_b \newline \implies &r = a - b = (\bar a + e_a) - (\bar b + e_b) = (\bar a - \bar b) + (e_a - e_b) \newline \implies & r = \bar r + (e_a - e_b) \newline \implies & \text{relative error of } r = \lvert \frac {e_a - e_b} {\bar r} \rvert \end{aligned} $$ 相比之前的 \(\text{relative error of } a \text{ or } b = \lvert \frac {e_a \text{ or } e_b} {a \text{ or } b} \rvert\)，分母减小了很多。
绝对误差¶
（计算同理）
数值稳定性¶
数值稳定性，简单来说，就是：
给入数据小误差，那么计算后的数据也应该小误差。
例子：\(I_n = \int_0^1 x^ne^x \mathrm x\)1¶
通过分部积分：
令 \(J_n = \frac 1 e I_n\)（归一化之后方便分析），则 $$ \begin{aligned} J_0 &= 1 - \frac 1 e\newline J_n &= 1 - n \frac 1 e I_{n-1} = 1 - n J_{n-1} \end{aligned} $$
Naive Approach¶
为了计算 \(J_{15}\)，我们需要不停地进行迭代，从而得出如此离谱的结果：
J(0): 0.6321206
J(1): 0.3678794
J(2): 0.2642412
J(3): 0.2072764
J(4): 0.1708944
J(5): 0.1455280
J(6): 0.1268320
J(7): 0.1121760
J(8): 0.1025920
J(9): 0.0766720
J(10): 0.2332800
J(11): -1.566080
J(12): 19.79296
J(13): -256.3085
J(14): 3589.319
J(15): -53838.78
这是因为，每一次计算中，我们实际上计算的是近似值：
（只考虑初始误差时，）这里的绝对误差是阶乘级（i.e. 甚至超过了指数级）增长的，相对误差更不必说。
Clever Approach¶
但是，还有另一种方法。我们可以反向计算：
我们容易证明：\(\frac {1}{n+1} = \int_0^1 x^ne^0 \dx \leq I_n = \int_0^1 x^ne^x \dx \leq \int_0^1 x^ne^1 \dx = \frac {e}{n+1}\)
因此：\(\frac{1}{e(n+1)} \leq J_n \leq \frac {1} {n+1}\)
并且：\(J_{n} = \frac 1 {n + 1} (1 - J_{n + 1})\)
我们可以先通过 \(J_{15}\) 的近似值，计算一下 \(J_1\) 试试看。
和准确值并无差别。
我们进而通过 \(J_{100}\) 的近似值，计算一下 \(J_0\) 到 \(J_{15}\)：
>>> max_val=100
>>> print_result(clever_J_upper)
clever_J_upper(0): 0.6321206
clever_J_upper(1): 0.3678794
clever_J_upper(2): 0.2642411
clever_J_upper(3): 0.2072766
clever_J_upper(4): 0.1708934
clever_J_upper(5): 0.1455329
clever_J_upper(6): 0.1268024
clever_J_upper(7): 0.1123835
clever_J_upper(8): 0.1009320
clever_J_upper(9): 0.09161229
clever_J_upper(10): 0.08387707
clever_J_upper(11): 0.07735223
clever_J_upper(12): 0.07177325
clever_J_upper(13): 0.0669477
clever_J_upper(14): 0.06273217
clever_J_upper(15): 0.05901754
>>> print_result(clever_J_lower)
clever_J_lower(0): 0.6321206
clever_J_lower(1): 0.3678794
clever_J_lower(2): 0.2642411
clever_J_lower(3): 0.2072766
clever_J_lower(4): 0.1708934
clever_J_lower(5): 0.1455329
clever_J_lower(6): 0.1268024
clever_J_lower(7): 0.1123835
clever_J_lower(8): 0.1009320
clever_J_lower(9): 0.09161229
clever_J_lower(10): 0.08387707
clever_J_lower(11): 0.07735223
clever_J_lower(12): 0.07177325
clever_J_lower(13): 0.0669477
clever_J_lower(14): 0.06273217
clever_J_lower(15): 0.05901754
可以发现 upper 和 lower 的值都在精度范围内相等。我们用了一个 \(J_{100}\) 的很粗糙的估计值，竟然算出了 \(J_0\) 到 \(J_{15}\) 的很精确的值！
这是因为，每一次计算中，我们都在减小误差：
可以发现，（只考虑初始误差时，）绝对误差从 max_val 到 n 是阶乘级下降的，相对误差更不必说。