误差

绪论

张瑞
中国科学技术大学数学科学学院

误差

误差

定义 1.
$\tilde x$为$x$的近似值，则称

\[E(x)=x-\tilde x \]

为误差，或绝对误差

\[R(x)=\frac{E(x)}{x}=\frac{x-\tilde x}{x} \]

为相对误差

有效位数

定义 2.
近似数的误差不超过某位的半个单位，则从这一位开始到第一个非0位的位数，称为有效位数

例 1. $\pi=3.1415926535\cdots$的近似值$3.14$， $3.141$, $3.14159$分别具有几位有效数字？

例 2. 若有数$x$，经过四舍五入后得到其近似值$\tilde x$，

\[\tilde x=\pm x_1.x_2\cdots x_n\times 10^m \]

其中$x_1\neq 0$，$x_1,\cdots,x_n$分别为$0,1,\cdots,9$中的一个数字，$m$为整数。则$\tilde x$的有效位数就是从最后一位到第一个非0位的位数，即$n$位有效数字。

例1：

  • 先看3.14中的最后一位。如果这是有效位，则真实值$\pi$应该落它的前后半个单位，即

    \[\pi\in[3.14-0.01\times\frac12, 3.14+0.01\times\frac12]=[3.135, 3.145] \]

    可以看到，这个条件是成立的。因此，3.14的有效位是从最后一位开始往前数到第一个非0位，即3位有效数字。

  • 对数3.141，如果最后一位是有效位，则应该有

    \[\pi\in[3.141-0.001\times\frac12,3.141+0.001\times\frac12]=[3.1405,3.1415] \]

    但$\pi>3.1415$，因此最后一位数1不是有效位。前面已经知道百分位是有效位， 因此，3.141同样只有3位有效数字。

例 3. 若$\tilde x$具有$n$位有效数字，估计其相对误差

解. 设

(1)
\[\tilde x=\pm x_1.x_2\cdots x_n\times 10^m \]

则$|x|>(x_1.x_2-0.1)\times 10^m$，因而

\[|R(x)|=|\frac{E(x)}{x}|\leq\frac{0.5\times 10^{-(n-1)}\times 10^m}{(x_1.x_2-0.1)\times 10^m} \]

简单点，可以用

\[|R(x)|=\leq\frac{0.5\times 10^{-(n-1)}\times 10^m}{x_1\times 10^m} \]

来估计相对误差

误差的来源

从实际问题到最后的计算结果产生的误差，主要有如下3类：

  • 原始误差，或称为模型误差，是在建模的过程中产生的误差
  • 截断误差，或称为方法误差，是数值方法产生的误差
  • 舍入误差，或称为计算误差，是由于计算机表达数据产生的误差

二进制表示

计算机中，用有限位数的二进制数来表达浮点数。 IEEE中，关于浮点数的标准

ieee-float

  • 由于机器字长的限制，计算机表达浮点数时，总会有偏差。
  • float类型的浮点数的精度范围是6~7位有效数字
  • double类型的精度范围是15~16位有效数字

例 4. 将十进制分数$\frac{11}{28}$，小数$0.4$表达为二进制数。将二进制数$110100.001011$表达为十进制数

解. 计算小数$0.4$，首先将小数部分一直乘2，积的整数部分顺序取出： 

0.4*2=0.8      取0      |                               
0.8*2=1.6      取1      |  顺
0.6*2=1.2      取1      |  序
0.2*2=0.4      取0      |  排
0.4*2=0.8      取0      |  列
0.8*2=1.6      取1      |
0.6*2=1.2      取1      |
0.2*2=0.4      取0      |

可以看出0110是循环，因此$0.4$的二进制是 

0.01100110……（循环0110）

十进制整数变为二进制数: 

11/2=5 余 1
5/2=2 余 1
2/2=1 余 0
1/2=0 余 1

则有 $11=(1011)_2$。类似有 $28=(11100)_2$。利用除法有 

            0.011 001 
        ———————— 
 11100 ) 1011.00 
          111 00 
         ————————— 
          100 000 
           11 100 
         ————————— 
              100 000 
               11 100 
             ————————— 
                  100

则有$\frac{11}{28}=0.011\dot00\dot1...$

对于二进制数$x_nx_{n-1}\cdots x_1x_0.y_1\cdots y_m$表达为十进制数为

\[x_n 2^n+x_{n-1}2^{n-1}+\cdots+x_1 2^1+x_0 2^0+y_1 2^{-1}+\cdots+y_m2^{-m} \]

例 5. 在单精度表达中，计算$10^6+10^{-6}$。最后的结果，会被表达成

\[1.000000000001\times 10^6 \]

但单精度只有小数点后6~7位的有效数字，也就是说，计算机只能表达成浮点数

\[1.0000000\times 10^6 \]

也就是说，

\[10^6+10^{-6}=10^6 \]

这种现象称为“大数吞小数”

误差的运算

若$\tilde x$, $\tilde y$为$x$, $y$的近似值。 用两个有误差的数($\tilde x$, $\tilde y$)做加减法后的值($\tilde x\pm\tilde y$) 与真实值($x$, $y$)做加减法得到的值($x+y$)之间的误差为(也就是误差的加减法)

\[(x\pm y)-(\tilde x \pm \tilde y)=(x-\tilde x) \pm (y-\tilde y) \]

有

\[E(x \pm y)=E(x) \pm E(y) \]

可以看到，加减后的绝对误差就是绝对误差的加减

相对误差

\[\begin{aligned} R(x + y)&=\frac{E(x)+ E(y)}{x+ y} =\frac{E(x)}{x}\frac{x}{x+ y}+\frac{E(y)}{y}\frac{y}{x+ y} \\ &=R(x)\frac{x}{x+ y}+R(y)\frac{y}{x+ y} \end{aligned} \]
  • 可以看到，当$x$与$y$同号时， $R(x+y)$在$R(x)$与$R(y)$之间。
  • 当$x$与$y$异号时，相对误差$R(x,y)$可能会变大。 特别当$x+y$接近0，或者说两个接近的数相减时，相对误差会变大，表现出来有效数字会减小。 如
    \[86.034-86.032=0.002 \]
    由5位有效数字，变为1位有效数字

乘法

\[xy-(\tilde x\tilde y)=x(y-\tilde y)+(x-\tilde x)\tilde y \]

有

\[E(xy)=xE(y)+E(x)\tilde y\approx xE(y)+yE(x) \]

相对误差

\[R(xy)=\frac{E(xy)}{xy}=\frac{E(y)}{y}+\frac{E(x)}{x}\frac{\tilde y}{y}\approx R(y)+R(x) \]

除法

\[\begin{aligned} \frac{x}{y}-\frac{\tilde x}{\tilde y} =&\frac{x\tilde y-\tilde x y}{y\tilde y}=\frac{x(\tilde y-y)+(x-\tilde x)y}{y\tilde y} \\ =&\frac{E(x)y-xE(y)}{y\tilde y} \end{aligned} \]

当$y$接近0时，绝对误差会被放大。

相对误差

\[R(\frac{x}{y})=\frac{E(\frac{x}{y})}{\frac{x}{y}}=\frac{E(x)}{x}\frac{y}{\tilde y}-\frac{E(y)}{\tilde y} \approx R(x)-R(y) \]

总结：

  1. 尽量避免除一个绝对值很小的数。如$h=0.01$，则编程时，优先使用$(\frac1h)^3$，而不是$\frac1{h^3}$。
  2. 尽量避免两个近似的数相减。 若$f(x_1)$与$f(x_2)$很接近时，应该换成其它的等价表达式，或者用近似公式
    \[f(x_1)-f(x_2)=f'(x_2)(x_1-x_2)+\frac{f''(x_2)}{2!}(x_1-x_2)^2+\cdots \]

例 6. 求$g(x)=10^7(1-\cos x)$在$2^{\circ}$处的近似值，用4位有效数字计算

解. 方法1： $\cos(2^{\circ})\approx 0.994$，则

\[g(2^{\circ})=10^7(1-\cos(2^{\circ})\approx 0.0006\times 10^7=0.6\times 10^4=6000 \]

（此时，数字有几位有效数字？）

方法2： $\cos(2^{\circ})=1-2\sin^2(1^{\circ})$，且$\sin(1^{\circ})\approx 0.01745$，则

\[\begin{aligned} g(2^{\circ})=10^7\times 2\times\sin^2(1^{\circ})\approx 0.0006090\times10^7 \\ =0.6090\times10^4=6090 \end{aligned} \]

（此时，数字有几位有效数字？）

真解：$0.60917298\times10^4$

例 7. 求$x^2-a x+1=0$的根，其中$a=2000000$， $a=200000000$

解. 根据公式

\[x_{1,2}=\frac{a\pm\sqrt{a^2-4}}2 \]

在双精度条件下，$a=2000000$时

\[x_2=5.00003807246685e-07 \]

此时，数据$x_2$的有效位数很少。（后面的3807246685是随机的）

若用$x_2=\frac1{x_1}$，则有

\[\tilde x_2=5.00000000000125e-07 \]

代入方程，可得

\[x_2^2-a x_2+1=-7.614493120033927e-06 \]
\[\tilde x_2^2-a \tilde x_2+1=0 \]

当$a=2000000$时

\[x_2=0.0 \]
\[\tilde x_2=5e-09 \]

作业

给出计算如下式子的方法，以达到相当的精度

  • $f(x)=(a+x)^n-a^n$，$x$接近$0$
  • $f(x)=\sin(a+x)-\sin(a)$，$x$接近$0$
  • $f(x)=x-\sqrt{x^2-a}$， $x>>a$

算法的稳定性

例 8. (舍入误差的作用) 求$\displaystyle I_n=\int_0^1\frac{x^n}{x+5}dx$

解. 由

\[I_{n}+5I_{n-1}=\int_0^1\frac{x^n+5x^{n-1}}{x+5}dx=\frac1n \]

可以构造2种算法:

\[\begin{aligned} I: &\quad I_n=\frac1n-5I_{n-1} \\ II: &\quad \tilde I_{n-1}=\frac15(\frac1n-\tilde I_n) \end{aligned} \]

可以用单精度，双精度，long double(C语言)来测试结果

运行结果表明：

  1. 算法I的计算会出现负值这种明显不合理的数值，而且随后数值开始振荡，而且振幅越来越大
  2. 算法I的振荡会随着表达精度（浮点数提高到双精度，再提高到long double）越高，出现的时间越晚
  3. 算法II不会出现算法I的不正常情况

两种算法均不存在方法上的误差，因此，不正常的现象只能是舍入误差引起的。

事实上，算法I中的数值的不正常现象的出现时间会随着表达精度的提高而变晚，就是一种典型的舍入误差导致的现象。

在算法I中，

  1. 若$I_0$带有误差$e_0$，则会放大5倍后，加入变量$I_1$，
  2. 然后又被放大5倍加入$I_2$. $...$
  3. 随着$n$越来越大，$e_0$也被放大的越来越大，最后起主导作用

而在算法II中，前一步的误差，缩小5倍后加入下一个值，因而最初的舍入误差的影响越来越小。

定义 3.
对舍入误差没有抑制作用的算法，称为不稳定的算法。

算法的稳定性与误差一样，是计算方法中非常重要的话题。

有时候，模型本身就是病态的

例 9. 解线性方程组

\[\begin{cases} x+ ay=1 \\ ax+y=0 \end{cases} \]

解. 可以看到，

  • 当$a=0.99$时，真解为$x=50.25$。
  • 当$a=0.991$时，真解为$x=55.81$。

$a$的变化只有大概$0.1\%$，但解的变化达到了$10\%$。解本身对输入参数很敏感。

如果一个模型的参数发生的很小的变化，会导致它的解产生很大的扰动，这样的模型称为病态的。

对于一个病态问题，

  • 需要更高精度浮点数的计算机程序来计算
  • 或者，想办法转换为等价、非病态的问题

程序作业

问题. 记级数(Hamming,1962)

\[\Phi(x)=\sum_{k=1}^\infty\frac{1}{k(k+x)} \]

编程计算$x$为$0,0.5,0.7,1$,$\sqrt2$,$10,100,300$时的值$\Phi(x)$， 并要求误差小于$10^{-6}$。


输出格式为 

x=0,       v=xxxxxxxxxxxx
x=0.5,     v=xxxxxxxxxxxx
x=0.7,     v=xxxxxxxxxxxx
x=sqrt(2), v=xxxxxxxxxxxx
x=10,      v=xxxxxxxxxxxx
x=100,     v=xxxxxxxxxxxx
x=300,     v=xxxxxxxxxxxx

计算量

例 10. 如何以最少的计算量来计算$x^{125}$

例 11. 如何以最少的计算量来计算

\[p_n(x)=a_0+a_1x+\cdots+a_nx^n \]

注意到

\[p_n(x)=((\cdots((a_nx+a_{n-1})x+a_{n-2})x+\cdots+a_2)x+a_1)x+a_0 \]

只要$n$次乘法，$n$次加法

目录

本节读完

例 12.

12.