---
url: /computational-physics/chapter-1/index.md
---
## Part 1 计算机表达的数

由于计算机都只有有限的内存，因此它只能够表达有限的数。这不仅仅包括实数，也包括整数。目前的经典计算机本质上都是通过记录二进制的数来记录所有的数的。也就是说，最基本的单元是所谓的**位 (bit)**，时髦的叫法是**比特**，它只有两个可能取值 $0$ 和 $1$。每 $8$ 个位构成一个**字节 (byte)**。

通常的整数由 $4$ 个 byte 构成，也就是说是 $32$ 位。当然，如果嫌这个不够，也可以用 $64$ 位 (也就是 $8$ 个 byte)。例如在 C 语言之中，`int` 型的整数就是 $4$ 个字节而所谓的 `long` 型的整数则是 $8$ 个字节。对任何一个 $4$ 字节 ($32\text{bit}$) 的非负整数都可以在二进制中表达为，

$$
n = \sum\_{i=0}^{31} b\_i 2^i
$$

正整数的范围因此从 $0$ 到 $2^{32}-1$。如果我们还希望表达负整数，我们只好牺牲一位来标记整数的正负，真正用于记录数字的位数减少到 $31$ 位。整数的最大特点是，只要在其定义域之内，机器对于整数的表达是严格的。而且整数之间的运算也是严格的 (只要结果不超出其定义域)。如果我们需要扩大整数的定义域，我们就需要启用更多位的整数。

对于实数，一般采用所谓的**浮点数 (floating point number)** 来表示。显然，由于实数构成不可数多的连续统，因此在计算机中能够表达的实数只能是其中的一个有限的部分。这就是浮点数系统。它保证了在任何一个实数的足够接近的邻域内总能找到一个浮点数来近似代表相应的实数。计算机中可以表达出来的那些数称为可表达的实数 (**representable real numbers**)，这个集合的大小取决于我们选择用多少位来表达一个实数。所有的那些不可表达的实数必须经过**舍入 (rounding)** 的操作，用它附近的一个可表达的实数来近似代替，相应产生的误差一般被称为**舍入误差**。特别地，我们可以定义一个在机器上能够表达的最小的正实数，它一般被称为**机器精度 (machine precision)**。它的定义可以取为：

$$
\epsilon\_M = \min {g \in A | 1 \oplus g > 1, g > 0}
$$

其中 $A$ 表示计算机中所有可以表示的实数而 $\oplus$ 表示计算机中实现的“加法” (里面包含了对结果的舍入)。显然，机器精度的数值依赖于我们在一个实数上愿意投入多少内存。例如，对于通常的单精度的实数 ($4$ 个字节) 的，$\epsilon\_M \sim 10^{-7}$ 而对于双精度的实数 ($8$ 个字节) 来说，$\epsilon\_M \sim 10^{-16}$。

上面的描述是对于浮点数系统的一个比较粗略的概述。如果说得更仔细一些的话，我们选取一个正的自然数 $\beta \in \mathbb{N}$ 为底来表达的话 (通常 $\beta=2$)，那么一般实数 $x$ 在以 $\beta$ 为底的系统中可以用下列整数部分和小数部分分别具有 $(n+1)$ 位和 $m$ 位的小数表达，

$$
x\_\beta = (-)^s \[x\_n x\_{n-1} \cdots x\_1 x\_0 . x\_{-1} x\_{-2} \cdots x\_{-m}] ;, \quad x\_n \ne 0
$$

其中 $s = 0,1$ 称为符号位，用以标记该数是正还是负。这个表达称为实数 $x$ (以 $\beta$ 为底) 的**位置表示 (positional representation)**。实数的位置表示的另外一个等价写法是

$$
x\_\beta = (-)^s \left\[ \sum\_{k=-m}^n x\_k \beta^k \right]
$$

事实上，更为经济与方便的计数方法是将实数提出一个公共的因子，

$$
x\_\beta = (-1)^s \cdot (0.a\_1 a\_2 \cdots a\_t) \cdot \beta^e = (-1)^s \cdot m \cdot \beta^{e-t} \tag{a}
$$

这个表示称为实数的**浮点数表示 (floating point representation)**，其中的 $s$ 仍然是符号位，整数 $m = a\_1 a\_2 \cdots a\_t$ 是一个 $t$ 位的整数称为**尾数 (mantissa)**，显然 $0 \le m \le \beta^t - 1$；整数 $e$ 则称为**指数 (exponent)**，它满足 $L \le e \le U$，一般选取 $L < 0$ 为一个负整数，$U > 0$ 为一个正整数。于是，对于一个占用 $N$ 位的实数而言，计算机中存储时会让符号位占一位，有效数字 $a\_i, i=1, \cdots, t$ 占 $t$ 位，指数 $e$ 则占据剩下的 $N - t - 1$ 位。

常见的情形是**单精度 (single precision)** 的浮点数 (占据 $32$ 位) 和**双精度 (double precision)** 的浮点数 (占据 $64$ 位)，它们除了符号位占一位之外，分别具有 $t=23$ 和 $t=52$ 位来存储有效数字，而存储指数的位数则分别是 $8$ 和 $11$ 位。一般来说数字 $0$ 会单独有一个特别的表示。因此，一个一般的浮点数系统 $\mathbb{F}$ 可以用它的基底 $\beta$、尾数的位数 $t$、指数 $e$ 的上下限 $L$ 和 $U$ 标记为，

$$
\mathbb{F}(\beta, t, L, U) = {0} \cup \left\[ x \in \mathbb{R} : x = (-)^s \beta^e \sum\_{k=1}^t a\_k \beta^{-k} \right]
$$

由于我们总是可以移动小数点并相应地改变指数 $e$ 的数值，因此为了保证数字表示的唯一性，一般假设 $a\_1 \ne 0$，否则我们可以向右移动小数点一位并减少 $e$ 一个单位即可。类似的，一般假设 $m \ge \beta^{t-1}$。经过这样约定的浮点数表示称为**正常化 (normalized)** 的浮点数表达。例如，在没有正常化的以 $10$ 为底的浮点数 $\mathbb{F}(10, 4, -1, 4)$ 系统中，实数 $1$ 可以有下列四种表达，

$$
0.1000 \cdot 10^1 ;, 0.0100 \cdot 10^2 ;, 0.0010 \cdot 10^3 ;, 0.0001 \cdot 10^4 ;,
$$

但是经过正常化后，只有第一个是符合要求的。容易证明，如果 $x \in \mathbb{F}$，那么我们必定有 $-x \in \mathbb{F}$。

容易发现，一个正常化表达的浮点数系统 $\mathbb{F}(\beta, t, L, U)$ 中的所有可表达的实数的绝对值一定介于 $x\_{\min}$ 和 $x\_{\max}$ 之间：

$$
x\_{\min} \equiv \beta^{L-1} \le |x| \le \beta^U (1 - \beta^{-t}) \equiv x\_{\max}
$$

为了能够表达更多 (特别是绝对值更小) 的数，一般会将上述正常化后的浮点数系统与 $e=L$ 时的放弃要求 $a\_1 \ne 0$ 的系统 (由于 $e=L$，因此这不会破坏数表示的唯一性) 合并。这可以产生在 $(-\beta^{L-1}, \beta^{L-1})$ 中的实数。这使得最小绝对值下降为 $\beta^{L-t}$。如果比这个还要小的数出现，机器会产生所谓的下溢 (underflow)。

浮点数的分布当然不是连续的。事实上也不是均匀的。假定有两个非零的浮点数 $x, y \in \mathbb{F}$，两者最近的距离一定介于 $\beta^{-1} \epsilon\_M |x|$ 和 $\epsilon\_M |x|$ 之间。其中的 $\epsilon\_M = \beta^{1-t}$ 就是前面引入的机器精度 (machine epsilon)，它实际上是最小的满足 $1 + \epsilon\_M > 1$ 的正实数。对于二进制来说，$\epsilon\_M \sim 10^{-7}$，如果取 $t=23$；$\epsilon\_M \sim 10^{-16}$，如果取 $t=52$。

如果我们考察相邻两个浮点数的相对误差，那么这是一个仅仅依赖于该数 mantissa 的函数，因此实际上是周期的。

$$
\frac{\Delta x}{x} = \frac{(-)^s \beta^{e-t}}{(-)^s m(x) \beta^{e-t}} = \frac{1}{m(x)}
$$

其中 $m(x)$ 表示实数 $x$ 的尾数部分 (一个整数)，由于 $\beta^{t-1} \le m(x) < \beta^t - 1$，因此上述的相对误差会在 $\[1/(\beta^t - 1), \beta^{1-t} = \epsilon\_M]$ 之间震荡。

对于那些不在系统 $\mathbb{F}$ 中的实数怎么办呢？计算机必须能够处理它们。最常见的操作称为**舍入 (rounding)**。这个操作是从一个一般的实数 $x \in \mathbb{R}$ 到某个浮点数系统 $\mathbb{F}$ 的映射，我们记为 $fl(\cdot)$。这个映射必须尽可能能够准确地反映原来的实数的大小，同时还不能改变实数的顺序。最简单的舍入方法就是利用大家熟悉的“四舍五入”原则。我们首先注意到任何的实数实际上都可以用公式 (a) 来表达出来，只要我们允许它的尾数可以有无穷多位。但是对一个给定的计算机浮点数系统，我们固定只能够保留 $t$ 位的尾数。但是我们总可以考察某个实数的尾数的下一位，$a\_{t+1}$ 并且进行四舍五入的操作：如果它比 $\beta/2$ 要小，就直接舍掉；如果它大或等于 $\beta/2$，就将前一位加 $1$：

$$
fl(x) = (-)^s (0.a\_1 a\_2 \cdots \tilde{a}*t) \cdot \beta^e ;, \quad \tilde{a}*t = \begin{cases} a\_t & a*{t+1} < \beta/2 \ a\_t + 1 & a*{t+1} \ge \beta/2 \end{cases}
$$

采用这样的处理之后，我们能够表达的实数范围大大扩展了。绝对值过大的实数仍然无法在计算机中表达。具体来说，如果 $|x| > x\_{\max}$，那么我们无法将其舍入到任何一个合理的可表达的实数。这时候计算机会产生一个**溢出 (overflow)**，程序的运行一般也会终止。类似的，如果某个实数的绝对值过小，$|x| < x\_{\min}$，这时候产生的一般称为**下溢 (underflow)**，此时往往系统会将其舍入为 $0$。

舍入造成的误差可以进行量化的估计。对于任意的 $x \in \mathbb{R}$，且满足 $x\_{\min} \le |x| \le x\_{\max}$，我们一定有，

$$
fl(x) = x(1 + \delta) ;, |\delta| \le \frac{1}{2} \beta^{1-t} = \frac{1}{2} \epsilon\_M
$$

因此我们有时候又称 $\epsilon\_M/2$ 为舍入误差单位 (**roundoff unit**)。这意味着一个实数与它的舍入值之间的误差是完全在控制之中的：

$$
E\_{rel}(x) \equiv \frac{|x - fl(x)|}{|x|} \le \frac{1}{2} \epsilon\_M ;, E(x) \equiv |x - fl(x)| \le \frac{1}{2} \beta^{e-t}
$$

从这个角度来说，浮点数的表示以及运算都是有误差的。浮点数比整数更为糟糕的是，即使是在其定义域内，浮点数也仅仅是某个实数的近似表示而不是严格表示，而且两个实数之间的运算的结果，即使它不超出定义域，通常也不是一个可表达的实数，需要将其经过舍入过程转换为一个可表达的实数。

一般来说，由于舍入问题的存在，计算机中的基本运算比如，加减法、乘法等等，并不满足原先数学中的结合律、分配律等基本规律。

例如，对于三个计算机中可表达的数 $a, b, c \in A$，在严格的数学中，$(a+b)+c = a+(b+c) = a+b+c$，但是在机器的运算中，$(a \oplus b) \oplus c$ 和 $a \oplus (b \oplus c)$ 一般是不相等的。

但是，舍入方法的实现应当尽量确保在一次基本运算后造成的舍入误差仅仅在一两个**机器精度**的范围之内。即使这一点得到满足，由于计算机在复杂的算法中往往需要进行 $N$ 多次的计算，因此我们必须考虑这种误差的传递和累计效应。

这使得我们在数值计算中对算法的稳定性的要求提高了。

## Part 2 算法稳定性

既然计算机中的实数都是利用不精确的浮点数来表达的，我们特别需要讨论算法的稳定性，或者说初始值中可能的误差，随着算法的运行是如何传递到最终的结果中去的。

我们首先来明确一下什么是算法。在计算数学中，从 $n$ 个已知的数据 $x \in \mathbb{R}^n$ 出发，最终获得最终 $m$ 个结果 $y \in \mathbb{R}^m$ 的一系列有限多次的操作之集合就构成了一个算法。为了获得一个具体问题的解可以有许多种算法。那么什么是一个好的算法呢。首先，一个好的算法首先必须是准确的 (当然，准确只是一个比较模糊的描述词语，具体的精度依赖于具体的问题)；也就是说，它可以给出问题的足够精确的近似解。其次，它必须是在计算资源方面比较快捷和有效的。当然正确性仍然是首要考虑。

既然每次计算机的操作往往都伴随着舍入误差，我们就需要分析这个误差是如何随着算法的发展而传递的。以一个简单的例子来说，假定我们需要计算三个实数 $a, b, c$ 的和，即 $a+b+c$，我们有两个算法来进行这个计算，即 $(a+b)+c$ 和 $a+(b+c)$。这两种算法对应的误差可以分析如下。

对 $(a+b)+c$ 而言，我们有，

$$
a \oplus b = (a+b)(1 + \varepsilon\_1) ;, (a \oplus b) \oplus c = \[(a+b)(1 + \varepsilon\_1) + c]\(1 + \varepsilon\_2) ;,
$$

我们可以忽略掉二阶的无穷小量得到这个算法的结果为，

$$
(a \oplus b) \oplus c = (a+b+c) \left\[ 1 + \frac{a+b}{a+b+c} \varepsilon\_2 + \varepsilon\_2 \right]
$$

显然，另外一种算法的误差估计就是将 $a$ 与 $c$ 互换，即，

$$
(b \oplus c) \oplus a = (a+b+c) \left\[ 1 + \frac{c+b}{a+b+c} \varepsilon\_2 + \varepsilon\_2 \right]
$$

因此，两种算法有一个部分是相同的，即 $\varepsilon\_2$ 的部分；但是另外一个部分的放大因子分别为 $(a+b)/(a+b+c)$ 和 $(b+c)/(a+b+c)$。注意，这两个因子可以相差很远！事实上，两者之比为，

$$
\left| \frac{a+b}{c+b} \right|
$$

这意味着，如果某两个之和恰好特别小，比如两个数几乎正负相消 (比如 $a+b \approx 0$)，那么先进行这两个数的加法的计算是更好的。事实上另外一种加法所产生的误差比首先进行几乎相消的两个数的加法要大很多倍。