---
url: /numerical-analysis/lesson-1-introduction/index.md
---
为什么要分析计算结果？

> 计算机算出的仅仅只是数值解，而不是精确解. 因此不能盲目采用，需要用不同的方法来具体分析.

误差类型：

* 模型误差：例如假设地球是一个球，忽略了一些次要因素.
* 测量误差
* 截断误差：用有限代替无限，用简单代替复杂导致的误差.
* 舍入误差：计算机能表示的数都是有限位的.

前两者并不是本课程预计要讨论的内容，但是我们要讨论它们产生的影响.

/Example/

> 用三种迭代算法计算 $3^{-n}$ ($n=1,2,\cdots,10$)：
>
> * $a\_0=0.99996$, $a\_n=a\_{n-1}/3$;
> * $b\_0=1$, $b\_1=0.33332$, $b\_n=4b\_{n-1}/3-b\_{n-2}/3$;
> * $c\_0=1$, $c\_1=0.33332$, $c\_n=10c\_{n-1}/3-c\_{n-2}/3$.
>
> ***
>
> 实际计算发现，后面两种有巨大的误差，给出的结果并不可靠.
>
> ```python
> a = [0.99996]
> b = [1.0, 0.33332]
> c = [1.0, 0.33332]
>
> for n in range(1, 11):
>     a.append(a[n - 1] / 3)
>
> for n in range(2, 11):
>     b.append(4 * b[n - 1] / 3 - b[n - 2] / 3)
>     c.append(10 * c[n - 1] / 3 - c[n - 2] / 3)
>
> print(f"{'n':>2} {'3^(-n)':>16} {'a_n':>16} {'b_n':>16} {'c_n':>16}")
> for n in range(1, 11):
>     print(f"{n:2d} {3.0 ** (-n):16.10e} {a[n]:16.10e} {b[n]:16.10e} {c[n]:16.10e}")
>
> ```
>
> 结果为：
>
> ```
> n           3^(-n)              a_n              b_n              c_n
>  1 3.3333333333e-01 3.3332000000e-01 3.3332000000e-01 3.3332000000e-01
>  2 1.1111111111e-01 1.1110666667e-01 1.1109333333e-01 7.7773333333e-01
>  3 3.7037037037e-02 3.7035555556e-02 3.7017777778e-02 2.4813377778e+00
>  4 1.2345679012e-02 1.2345185185e-02 1.2325925926e-02 8.0118814815e+00
>  5 4.1152263374e-03 4.1150617284e-03 4.0953086420e-03 2.5879159012e+01
>  6 1.3717421125e-03 1.3716872428e-03 1.3517695473e-03 8.3593236214e+01
>  7 4.5724737083e-04 4.5722908093e-04 4.3725651578e-04 2.7001773438e+02
>  8 1.5241579028e-04 1.5240969364e-04 1.3241883859e-04 8.7219470251e+02
>  9 5.0805263425e-05 5.0803231215e-05 3.0806279531e-05 2.8173097636e+03
> 10 1.6935087808e-05 1.6934410405e-05 -3.0645734898e-06 9.1003009778e+03
> ```
>
> 后两列很明显不正常.

原理在于，计算机中 $\beta$ 进制非零浮点数尾数的一般形式为
$$
\pm0.d\_1d\_2\cdots d\_t\times \beta^J
$$
其中，尾数 $0.d\_1d\_2\cdots d\_t$ 中 $0\leqslant d\_t\leqslant \beta-1$ ($k=1,2,\cdots,t$)，$d\_t>0$. $t$ 是尾数的位数，称为「字长」；整数 $J$ 称为「阶」，满足 $L\leqslant J\leqslant U$，$L,U$ 有很多种标准值.

对于二进制，给定一套 $(t,L,U)$，就能给出一个有限的数集 (包含零)，总共数量有：
$$
2^{t-1}(U-L+1)+1
$$
这个数集中的非零数，对称地、==不等距地== 分布在区间 $\[-M,-m]$ 和 $\[m,M]$ 上，其中
$$
m=2^{t-1},\quad M=2^U(1-2^{-t})
$$
分别是系统所能表示的最小和最大正浮点数. 对于 $j=L, L+1,\cdots,U$，区间 $\[2^{t-1},2^j)$ 内部等距分布 $2^{t-1}$ 个数；但是区间的长度是会变化的，因此在整个实数轴上，浮点数的分布并不等距.

对于实数 $x$，相应的浮点数记为 $\text{fl}(x)$：

* 若 $x=0$，则 $\text{fl}(x)=0$;
* 若 $m\leqslant |x|\leqslant M$，则 $\text{fl}(x)$ 为整个浮点数集 $F$ 中最接近 $x$ 的那个浮点数;
* 若不符合上述条件，那么就会「上溢」或者「下溢」.

/Theorem/

> 设实数 $x$ 满足 $m\leqslant |x|\leqslant M$，则存在实数 $\delta$ 满足 $|\delta|\leqslant 2^{-t}$，使得
> $$
> \text{fl}(x)=x(1+\delta)
> $$
> ($\delta$ 的上界称为机器精度 $\varepsilon\_{\text{mach}}$.)

有了这个规则之后，可以建立起浮点数的「加减乘除」.

* 加减法先对阶，后运算，再舍入.
* 乘除法先运算，再舍入.
* 不在浮点数系中的数做舍入 (最靠近) 处理.

因此在计算机系统中加法结合律实际上并不成立，当然这个误差是很小的.

误差的传播和估计：对于多元 (可微) 函数 $f=f(x\_1,x\_2,\cdots,x\_d)$，设自变量 $X\equiv(x\_1,x\_2,\cdots,x\_d)$ 的近似值是 $X\_A$，那么函数值的误差为
$$
f(X)-f(X\_A)=\sum\_{j=1}^d(x\_j-x\_{jA})\int\_0^1\frac{\partial f(\theta X-(1-\theta)X\_A)}{\partial x\_j}\mathrm{d}\theta
$$
误差估计：
$$
|f(X)-f(X\_A)|=\sum\_{j=1}^d(x\_j-x\_{jA})\int\_0^1\left|\frac{\partial f(\theta X-(1-\theta)X\_A)}{\partial x\_j}\right|\mathrm{d}\theta
$$
有效数字：一个实数 $x$ 的近似值 $x\_A$ 用十进制表示为
$$
x\_A=\pm 10^k\cdot0.d\_1d\_2\cdots d\_j
$$
如果
$$
|x-x\_A|\leqslant 0.5\times 10^{k-j}
$$
那么 $x\_A$ 为 $x$ 的具有 $j$ 位的 (十进制) 有效数字近似值.

/Definition/ (条件数)

> 称 $C=C(x)$ 为 $f$ 在点 $x$ 的条件数，若：
> $$
> \frac{|f(x)-f(x\_A)|}{|f(x)|}\leqslant C(x)\frac{|x-x\_A|}{|x|}
> $$
> 反映函数随自变量的变化程度，是函数的本身性质；如果函数本身可微，那么条件数就接近 $|xf'(x)|/|f(x)|$.
>
> 条件数大的问题称为病态问题；反之为良态问题.

如果一个算法可以保证初始小误差不造成结果的重大影响，那么这是一个稳定的算法 —— 在实际应用中我们只应该使用稳定的算法.

误差的避免：

* 防止接近零的数作为除数，放大误差

* 防止大数「吃掉」小数，加减法运算中的大数的误差本身就会掩盖掉小数

* 防止相近的两个数相减，损失有效数字，一个例子是分子的有理化，比如
  $$
  8-\sqrt{63}=\frac{1}{8+\sqrt{63}}
  $$
  有效数字从 1 位变成 2 位.

* 减少运算次数

  计算量是另一个需要考量的指标，我们一般把一个算法所需要的乘除运算总次数称为计算量，单位是 flop (也就是单次乘除法计算，floating point operation).

  比如秦九韶算法：
  $$
  a\_5x^5+a\_4x^4+a\_3x^3+a\_2x^2+a\_1x^1+a\_0=((((a\_5x+a\_4)x+a\_3)x+a\_2)x+a\_1)x+a\_0
  $$
  从 15 次 flop 减少到 5 次.

除此之外还应该考虑存储量，即算法占用多大内存. 比如要存储前两步的结果，或者只用存储上一步的结果.