19.8.2 计算机计算中的数值问题

19.8.2.1 引言、误差类型

计算机计算的一般性质与手工计算基本相同, 然而某些方面需要特别注意, 因为精度来自数的表示和关于计算机误差的判断. 更进一步, 计算机要比人类手工能做的实施多得多的计算步骤.

因此, 存在如何影响和控制误差的问题, 例如在数学上等价的方法中选用最适当的数值方法.

在后面的讨论中用到如下符号,其中 $x$ 表示在大多情况下未知的量的准确值, $\tilde{x}$ 表示 $x$ 的近似值:

绝对误差 $\left|\mathrm{\Delta }x\right|=\left|x-\tilde{x}\right|$ .(19.268)

相对误差 $\left|\frac{\mathrm{\Delta }x}{x}\right|=\left|\frac{x-\tilde{x}}{x}\right|$ .(19.269)

记号

$$\begin{array}{}\text{(19.270)}& \epsilon \left(x\right)=x-\tilde{x}\text{ 和 }{\epsilon }_{\text{rel }}\left(x\right)=\frac{x-\tilde{x}}{x}\end{array}$$

也经常用到.

19.8.2.2 规范化十进制数和舍入

1. 规范化十进制数

每个实数 $x\ne 0$ 可表示为形如

$$\begin{array}{}\text{(19.271)}& x=\pm 0.b_1{b}_2\cdots {10}^E\left(b_1\ne 0\right)\end{array}$$

的十进制数. 这里由数字 $b_i\in \{0,1,2,\cdots ,9\}$ 构成的 $0.b_1{b}_2\cdots$ 称为尾数. 数 $E$ 为整数,是所谓关于基 10 的指数. 因为 $b_1\ne 0,\left(19.271\right)$ 称为正规十进制数.

因为真实的计算机只能处理有限多的字节,故必须限制尾数数字的固定数目 $t$ 和指数 $E$ 的固定范围. 故形如 (19.271) 的数根据舍入 (在实际计算中常用) 得到

$$\begin{array}{}\text{(19.272)}& \tilde{x}=\{\begin{array}{ll}\pm 0\cdot b_1{b}_2\cdots b_t\cdot {10}^E,& b_{t+1}<5\left(\text{ 舍 }\right),\\ \pm \left(0\cdot b_1{b}_2\cdots b_t+{10}^{-t}\right){10}^E,& b_{t+1}\ge 5\left(\lambda \right),\end{array}\end{array}$$

由舍入引起的绝对误差为

$$\begin{array}{}\text{(19.273)}& \left|\mathrm{\Delta }x\right|=\left|x-\tilde{x}\right|\le 0.5\cdot {10}^{-t}{10}^E.\end{array}$$

2. 基本运算和数值计算

每个数值过程都是一系列基本运算. 特别用有限位浮点表示提出问题. 这里给出简要综述. 设 $x$ 和 $y$ 是非零的同号规范化无误差浮点数:

$$\begin{array}{}\text{(19.274a)}& x=m_1{B}^{E_1},y=m_2{B}^{E_2}\end{array}$$$$\begin{array}{}\text{(19.274b)}& m_i=\sum _{k=1}^t{a}_{-k}^{\left(i\right)}{B}^{-k},a_{-1}^{\left(i\right)}\ne 0,\end{array}$$$$\begin{array}{}\text{(19.274c)}& a_{-k}^{\left(i\right)}=0\text{或 1 或}\cdots \text{或}B-1,k>1\left(i=1,2\right)\text{.}\end{array}$$

(1) 加法 若 $E_1>E_2$ ,因为正规化仅允许左移,则公共指数变为 $E_1$ . 随后尾数相加.

若

$$\begin{array}{}\text{(19.275a)}& B^{-1}\le \left|m_1+m_2{B}^{-\left(E_1-E_2\right)}\right|<2\end{array}$$

和

$$\begin{array}{}\text{(19.275b)}& \left|m_1+m_2{B}^{-\left(E_1-E_2\right)}\right|\ge 1,\end{array}$$

则将十进制小数点向左移一位而指数增加 1.

$0.9604\cdot {10}^3+0.5873\cdot {10}^2=0.9604\cdot {10}^3+0.05873\cdot {10}^3=1.01913\cdot {10}^3=0.1019\cdot {10}^4.$

(2) 减法 如同在加法的情况均衡指数, 随后尾数相减. 若

$$\begin{array}{}\text{(19.276a)}& \left|m_1-m_2{B}^{-\left(E_1-E^2\right)}\right|<1-B^{-t},\end{array}$$

以及

$$\begin{array}{}\text{(19.276b)}& \left|m_1-m_2{B}^{-\left(E_1-E^2\right)}\right|<B^{-t},\end{array}$$

则将十进制小数点右移 $t$ 的最大值位,而指数相应减少.

$◼$ $0.1004\cdot {10}^3-0.9988\cdot {10}^2=0.1004\cdot {10}^3-0.09988\cdot {10}^3=0.00052\cdot {10}^3=0.5200\cdot {10}^0.$

此例显示了减法的临界情况. 因为位数有限 (这里是 4), 从右边引进零代替准确字符.

(3) 乘法 指数相加而尾数相乘, 若

$$\begin{array}{}\text{(19.277)}& m_1{m}_2<B^{-1},\end{array}$$

则十进制小数点向右移一位, 且指数减少 1.

$\left(0.1004\cdot {10}^3\right)\cdot \left(0.2504\cdot {10}^5\right)=0.07952704\cdot {10}^8=0.7953\cdot {10}^7.$

(4)除法 指数相减而尾数相除. 若

$$\begin{array}{}\text{(19.278)}& \frac{m_1}{m_2}\ge 1\end{array}$$

则十进制小数点向左移一位, 且指数增加 1.

$\left(0.3176\cdot {10}^3\right)/\left(0.2504\cdot {10}^5\right)=1.2683706\cdot {10}^{-2}=0.1268\cdot {10}^{-1}.$

(5) 结果的误差 在假定无误差项的四种基本运算中, 结果的误差是舍入误差. 对于位置为 $t$ 基为 $B$ 的数,相对误差的上限为

$$\begin{array}{}\text{(19.279)}& \frac{B}{2}{B}^{-t}\text{.}\end{array}$$

(6) 减法相消 如上所述, 几乎相等的浮点数的减法是临界运算. 若有可能, 应通过改变运算阶或利用某种等式来避免这种情况.

$$x=\sqrt{1985}-\sqrt{1984}=0.4455\cdot {10}^2-0.4454\cdot {10}^2=0.1\cdot {10}^{-1}\text{或}$$$$x=\sqrt{1985}-\sqrt{1984}=\frac{1985-1984}{\sqrt{1985}+\sqrt{1984}}=0.1122\cdot {10}^{-1}.$$

19.8.2.3 数值计算的精度

1. 误差类型

数值方法有误差. 有几类误差, 最后结果的总误差正是由这些误差积累的 (图 19.18).

2. 输入误差

(1)输入误差的概念 输入误差是由不准确的输入数据产生的误差. 输入数据的轻微不准确称为扰动. 确定输入数据误差称为误差计算的直接问题. 其反问题如下: 输入数据可有多大的误差能够保证最终的输入误差不超过可以接受的允许值. 在相当复杂的问题里估计输入误差是非常困难且通常几乎是不可能的. 一般对实值函数 $y=f\left(\underset{―}{\mathbf{x}}\right),\underset{―}{\mathbf{x}}={(x_1,x_2,\cdots ,x_n)}^{\mathrm{T}}$ ,若对 $y=f\left(\underset{―}{\mathbf{x}}\right)=f\left(x_1,x_2,\cdots ,x_n\right)$ 应用带线性余项的泰勒公式 (参见第 630 页 7.3.3.3), ${\xi }_1,{\xi }_2,\cdots ,{\xi }_n$ 表示中间值, ${\tilde{x}}_1,{\tilde{x}}_2,\cdots ,{\tilde{x}}_n$ 表示 $x_1,x_2,\cdots ,x_n$ 的近似值,则输入误差的绝对值为

$$\left|\mathrm{\Delta }y\right|=\left|f\left(x_1,x_2,\cdots ,x_n\right)-f\left({\tilde{x}}_1,{\tilde{x}}_2,\cdots ,{\tilde{x}}_n\right)\right|$$$$=\left|\sum _{i=1}^n\frac{\partial f}{\partial x_i}\left({\xi }_1,{\xi }_2,\cdots ,{\xi }_n\right)\left(x_i-\widetilde{x_i}\right)\right|$$$$\begin{array}{}\text{(19.280)}& \le \sum _{i=1}^n\left(\underset{x}{max}\left|\frac{\partial f}{\partial x_i}\left(\underset{―}{\mathbf{x}}\right)\right|\right)\left|\mathrm{\Delta }{x}_i\right|,\end{array}$$

近似值是扰动了的输入数据. 这里也考虑高斯误差传播定律 (参见第 1114 页 16.4.2.1).

(2)简单算术运算的输入误差 已知简单算术运算的输入误差. 对四种基本运算用 (19.268) (19.270) 的记号:

$$\begin{array}{}\text{(19.281)}& \epsilon \left(x\pm y\right)=\epsilon \left(x\right)\pm \epsilon \left(y\right),\end{array}$$$$\begin{array}{}\text{(19.282)}& \epsilon \left(xy\right)=y\epsilon \left(x\right)+x\epsilon \left(y\right)+\epsilon \left(x\right)\epsilon \left(y\right),\end{array}$$$$\begin{array}{}\text{(19.283)}& \epsilon \left(\frac{x}{y}\right)=\frac{1}{y}\epsilon \left(x\right)-\frac{x}{y^2}\epsilon \left(y\right)+\epsilon \text{的高阶项,}\end{array}$$$$\begin{array}{}\text{(19.284)}& {\epsilon }_{\text{rel }}\left(x\pm y\right)=\frac{{\epsilon }_{\text{rel }}\left(x\right)\pm {\epsilon }_{\text{rel }}\left(y\right)}{x\pm y},\end{array}$$$$\begin{array}{}\text{(19.285)}& {\epsilon }_{\text{rel }}\left(xy\right)={\epsilon }_{\text{rel }}\left(x\right)+{\epsilon }_{\text{rel }}\left(y\right)+{\epsilon }_{\text{rel }}\left(x\right){\epsilon }_{\text{rel }}\left(y\right),\end{array}$$$$\begin{array}{}\text{(19.286)}& {\epsilon }_{\text{rel }}\left(\frac{x}{y}\right)={\epsilon }_{\text{rel }}\left(x\right)+{\epsilon }_{\text{rel }}\left(y\right)+\epsilon \text{的高阶项.}\end{array}$$

公式表明: 对于乘法和除法, 输入数据的相对误差小, 导致结果的相对误差也小. 对于加法和减法,若 $\left|x\pm y\right|\ll \left|x\right|+\left|y\right|$ ,相对误差可能非常大.

3. 方法的误差

(1)方法误差的记号 方法误差源于理论上连续的现象作为极限以不同的方式被数值逼近的事实. 因此, 在极限过程中有截断误差 (例如在迭代法中) 及在用有限离散系 (例如数值积分) 逼近连续现象时的离散误差. 方法误差与输入和舍入误差无关, 因此, 仅在关系到应用解法的方法论时研究方法误差.

(2) 应用迭代法 若使用迭代法, 可能出现两种情况: 得到问题的正确解或错误解. 也可能尽管有解但不能用迭代法得到.

为使迭代法更清晰安全, 应考虑如下建议:

a) 为避免 “无穷迭代”, 若步数超过预定值即停止过程 (即尚未达到要求的精度便停止).

b) 应在屏幕上以数值或者图表的形式跟踪中间结果的位置.

c) 应该用到解的所有已知性质如梯度、单调性等.

d) 应研究变量和函数计量的可能性.

e) 应通过改变步长、截断条件、初始值等进行多种试验.

4. 舍入误差

产生舍入误差是因为中间结果被舍入. 这对按精度要求判断数学方法时有本质的重要性. 舍入误差与输入误差和方法误差一起决定给定的方法是强稳定、弱稳定或不稳定. 若总误差随着步数增加分别减少、有相同的阶或增加, 便发生强稳定、弱稳定或不稳定.

在不稳定性方面, 我们区别舍入误差和离散误差 (数值不稳定) 以及理论上准确的计算中初始数据误差 (自然不稳定) 的灵敏度. 若数值不稳定不大于自然不稳定, 则计算过程是合适的.

对于舍入误差的局部误差传播, 即从一个计算步到下一步的误差传递, 可使用在输入误差中用过的同样的估计过程.

5. 数值计算的例子

上述某些问题用数值例子来说明.

$◼\mathbf{A}$: 二次方程的根 带实系数 $a,b,c$ 的二次方程 $a{x}^2+bx+c=0,D=$ $b^2-4ac\ge 0$ (实根). 临界状态为

**a) $|\ast \ast 4ac|\ll b^2$ 和 b) $4ac\approx b^2$ . 推荐程序:

**i) $x_1\ast \ast =-\frac{b+\mathrm{sign}\left(b\right)\sqrt{D}}{2a},x_2=\frac{c}{a{x}_1}$ (韦达根定理,参见第 56 页 1.6.3.1,3.).

ii) 用直接法难免把 $D$ 化零. 因为 $\left|b\right|\gg \sqrt{D}$ 成立,将发生减法抵消,除非 $\left(b+\mathrm{sign}\left(b\sqrt{D}\right)\right)$ 中误差不是太大.

$◼\mathbf{B}$: $h\le r$ 的薄锥壳的体积 因为 $\left(r+h\right)\approx r,V=4\pi \frac{{(r+h)}^3-r^3}{3}$ 存在减法消去的情况. 而在等式 $V=4\pi \frac{3r^{2}h+3r{h}^2+h^3}{3}$ 中则没有这个问题.

$◼\mathbf{C}$: 求和 $S=\sum _{k=1}^{\mathrm{\infty }}\frac{1}{k^2+1}\left(S=1.07667\cdots \right)$ 要求有三位有效数字的精度. 用 8 位数字进行计算,大约需要加 6000 项. 在作恒等变换 $\frac{1}{k^2+1}=\frac{1}{k^2}-\frac{1}{k^2\left(k^2+1\right)}$ 后, 成立

$$S=\sum _{k=1}^{\mathrm{\infty }}\frac{1}{k^2}-\sum _{k=1}^{\mathrm{\infty }}\frac{1}{k^2\left(k^2+1\right)}\text{ 及 }S=\frac{{\pi }^2}{6}-\sum _{k=1}^{\mathrm{\infty }}\frac{1}{k^2\left(k^2+1\right)}.$$

通过这一变换后, 则只需考虑 8 项.

$◼\mathbf{D}$: 避免 $\frac{0}{0}$ 的状态 当 $x=y=0$ 时,函数 $z=\left(1-\sqrt{1+x^2+y^2}\right)\frac{x^2-y^2}{x^2+y^2}$ . 分子和分母同时乘以 $\left(1+\sqrt{1+x^2+y^2}\right)$ 即可避免这一状态.

$◼\mathbf{E}$: 不稳定递推过程的例子 若满足条件 $\left|\frac{a}{2}\pm \sqrt{\frac{a^2}{4}+b}\right|<1$ ,则一般形式的算法 $y_{n+1}=a{y}_n+b{y}_{n-1}\left(n=1,2,\cdots \right)$ 是稳定的. 特殊情况 $y_{n+1}=-3y_n+$ $4y_{n-1}\left(n=1,2,\cdots \right)$ 是不稳定的. 若 $y_0$ 和 $y_1$ 有误差 $\epsilon$ 和 $-\epsilon$ ,则对 $y_2,y_3,y_4,y_5$ , $y_6,\cdots$ 误差为 $7\epsilon ,-25\epsilon ,103\epsilon ,-409\epsilon ,1639\epsilon ,\cdots$ ,该过程对于参数 $a=-3$ 与 $b=4$ 是不稳定的.

$◼\mathbf{F}$: 微分方程的数值求积 数值求解一阶常微分方程

$$\begin{array}{}\text{(19.287)}& y^{\prime }=f\left(x,y\right),\text{ 其中 }f\left(x,y\right)=ay,\end{array}$$

其初值用 $y\left(x_0\right)=y_0$ 表示.

a) 天然不稳定 准确解 $y\left(x\right)$ 有准确初值 $y\left(x_0\right)=y_0$ ,设 $u\left(x\right)$ 为扰动初值的解. 不失一般性, 设扰动解形如:

$$\begin{array}{}\text{(19.288a)}& u\left(x\right)=y\left(x\right)+\epsilon \eta \left(x\right),\end{array}$$

其中 $\epsilon \left(0<\epsilon <1\right)$ 为参数,而 $\eta \left(x\right)$ 是所谓扰动函数. 考虑 $u^{\prime }\left(x\right)=f\left(x,u\right)$ 从泰勒展开式得到 (参见第 630 页 7.3.3.3)

$$\begin{array}{}\text{(19.288b)}& u^{\prime }\left(x\right)=f\left(x,y\left(x\right)+\epsilon \eta \left(x\right)\right)=f\left(x,y\right)+\epsilon \eta \left(x\right)f_y\left(x,y\right)+\text{ 高阶项. }\end{array}$$

这意味着微分变差方程

$$\begin{array}{}\text{(19.288c)}& {\eta }^{\prime }\left(x\right)=f_y\left(x,y\right)\eta \left(x\right).\end{array}$$

由 $f\left(x,y\right)=ay$ ,问题的解为

$$\begin{array}{}\text{(19.288d)}& \eta \left(x\right)={\eta }_0{\mathrm{e}}^{a\left(x-x_0\right)},\text{ 其中 }{\eta }_0=\eta \left(x_0\right).\end{array}$$

对 $a>0$ ,即便是小的初始扰动 ${\eta }_0$ 也导致无限增长的扰动 $\eta \left(x\right)$ . 故为天然不稳定.

b) 梯形公式的误差研究 当 $a=-1$ 时,稳定的微分方程 $y^{\prime }\left(x\right)=-y\left(x\right)$ 有准确解

$$\begin{array}{}\text{(19.289a)}& y\left(x\right)=y_0{\mathrm{e}}^{a\left(x-x_0\right)},\text{ 其中 }{y}_0=y\left(x_0\right).\end{array}$$

梯形公式为

$$\begin{array}{}\text{(19.289b)}& {\int }_{x_i}^{x_{i+1}}y\left(x\right)\mathrm{d}x\approx \frac{y_i+y_{i+1}}{2}h,\text{ 其中 }h=x_{i+1}-x_i.\end{array}$$

对给定微分方程用上述公式, 成立

$${\tilde{y}}_{i+1}={\tilde{y}}_i+{\int }_{x_i}^{x_{i+1}}\left(-y\right)\mathrm{d}x={\tilde{y}}_i-\frac{{\tilde{y}}_i+{\tilde{y}}_{i+1}}{2}h,{\tilde{y}}_{i+1}=\frac{2-h}{2+h}{\tilde{y}}_i,$$$$\begin{array}{}\text{(19.289c)}& {\tilde{y}}_i={\left(\frac{2-h}{2+h}\right)}^i{\tilde{y}}_0\end{array}$$

其中 $x_i=x_0+ih$ ,即对 $0\le h<2$ 有 $i=\left(x_i-x_0\right)/h$ ,得到

$${\tilde{y}}_i={\left(\frac{2-h}{2+h}\right)}^{\left(x_i-x_0\right)/h}{\tilde{y}}_0={\tilde{y}}_0{\mathrm{e}}^{c\left(h\right)\left(x_i-x_0\right)},$$$$\begin{array}{}\text{(19.289d)}& c\left(h\right)=\frac{\ln \left(\frac{2-h}{2+h}\right)}{h}=-1-\frac{h^2}{12}-\frac{h^4}{80}-\cdots .\end{array}$$

若 ${\tilde{y}}_0=y_0$ ,则 ${\tilde{y}}_i<y_i$ ,且对 $h\to 0,{\tilde{y}}_i$ 也趋向于准确解 $y_0{\mathrm{e}}^{-\left(x_i-x_0\right)}$ .

c) 在 b) 中的输入误差 设准确和近似的初值相同. 现研究当 ${\tilde{y}}_0\ne y_0$ 及 $\left|{\tilde{y}}_0-y_0\right|<{\epsilon }_0$ 时的性态.

因为 $\left({\tilde{y}}_{i+1}-y_{i+1}\right)\le \frac{2-h}{2+h}\left({\tilde{y}}_i-y_i\right)$ ,有

$$\begin{array}{}\text{(19.290a)}& \left({\tilde{y}}_{i+1}-y_{i+1}\right)\le {\left(\frac{2-h}{2+h}\right)}^{i+1}\left({\tilde{y}}_0-y_0\right),\end{array}$$

故 ${\epsilon }_{i+1}$ 最多和 ${\epsilon }_0$ 同阶,且该方法关于初值是稳定的. 应该提到,在用辛普森方法求解上述微分方程时引进了人为的不稳定. 此时,对 $h\to 0$ ,得到通解如:

$$\begin{array}{}\text{(19.290b)}& {\tilde{y}}_i=C_1{\mathrm{e}}^{-x_i}+C_2{(-1)}^i{\mathrm{e}}^{x_i/3}.\end{array}$$

问题是该数值解法使用了比相应的微分方程的阶更高阶的差分.