Skip to main content Link Search Menu Expand Document (external link)
Fiddie
目录

上机作业基本要求

问题要回答清楚、结合上机试验得到的结果以及现成的理论进行解释．

答疑：“精度丢失定理”如何理解？

“精度丢失定理”如何理解？

先体会一下精度丢失

在计算机中，运行如下 python 代码：

x = 0.1
y = 0.4
print(y-x)

得到的结果是0.30000000000000004，而不是0.3！

注： 二进制转十进制的 python 代码：

def bin_to_dec(b):
    p=b.index('.')
    return int(b[:p],2)+int(b[p+1:],2)/2**len(b[p+1:])

原因？

这是因为，在IEEE双精度浮点数系统中，没有办法精确存储 0.3，只能存储它的一个舍入近似值．

将0.4和0.1的二进制形式按实际展开，末尾补零相减，结果如下：

 0.01100110011001100110011001100110011001100110011001101000
-0.00011001100110011001100110011001100110011001100110011010
=0.01001100110011001100110011001100110011001100110011001110

因为双精度浮点数的机器精度是

\[\varepsilon_{\mathrm{Mach}}=2^{-52} \approx 2.2\times 10^{-16},\]

所以保留 52 位有效数字之后会得到下面的规格化浮点数：

\[0.1001100110011001100110011001100110011001100110011010 \times 10^{-2}.\]

把上面的数转成十进制，就是0.30000000000000004．然而真实值写成二进制应该是 $0.01001100\cdots$（1100循环）．

因为 $1-\dfrac{x}{y} = 0.75 > 2^{-1}$，根据精度丢失定理，丢失掉 1 个有效的二进制位．

如何判断精度丢失 1 位？

假设我们知道真实值 $x$ 的规格化浮点数，以及按照算法得出来的数字 $\overline{x}$ 的规格化浮点数．为了在计算机中验证一个算法的精度丢失不超过1位，只需要判断相对误差是否满足下面的不等式：

\[\left\vert\dfrac{\overline{x}-x}{x}\right\vert \le \varepsilon_{\mathrm{Mach}}.\]

就比如，采用前面的例子，$x = 0.3$，$\overline{x}=0.30000000000000004$，于是

\[\left\vert\dfrac{\overline{x}-x}{x}\right\vert = \dfrac{4\times 10^{-17}}{0.3} \approx 1.33 \times 10^{-16} < \varepsilon_{\mathrm{Mach}}.\]

所以在计算机中计算 $0.1+0.2$ 会让精度丢失不超过 1 位（其实就恰好丢失了 1 位）．

第1章上机作业

问题1： 设计算法找出你所使用的计算机的机器精度、下溢值和上溢值，编程验证所设计算法的正确性．

问题2： 设计算法计算 $y= x−\sin x$，使得有效位的丢失最多 1 位．编程验证所设计算法的正确性．

问题3： 计算积分 $\displaystyle \int_0^1x^n\mathrm{e}^x\mathrm{d}x(n\ge 0)$，可用分部积分法计算得 $y_{n+1}=\mathrm{e}-(n+1)y_n$，编程验证按上述公式设计的算法不稳定．

问题4： 考虑由 $x_0=1$，$x_1=\dfrac{1}{3}$，$x_{n+1}=\dfrac{13}{3}x_n-\dfrac{4}{3}x_{n-1}$ 归纳定义的实数序列实现的算法是否稳定？如果将初值改为 $x_0=1$，$x_1=4$，数值稳定吗？

作业评析

本次作业评分标准为 $x+y+z$ ，其中：

1．基准分 $x$：
① 上机报告的基本要素都有，基准分为 10 分．
② 上机报告缺了任何一个部分，基准分为 9 分．
③ 上机报告缺了至少两个部分，基准分为 7 分．
④ 雷同或抄袭：基准分为 0 分．

2．减分：

减分点个数得分 $y$
$0\sim 1$$0$
$2$$-1$
$3$$-2$
$\ge 4$$-3$

3．加分 $z$：
如果有亮点（详见下面的各个问题），加一定的分数．

问题1

我们采用IEEE双精度浮点数(double类型)．

1．机器精度(machine epsilon) 是满足

\[fl(1+\varepsilon) > 1\]

的最小正机器数 $\varepsilon$．根据浮点数存储方式，可以知道机器精度是

\[2^{-52}\approx 2.2\times 10^{-16}.\]

2．下溢值(underflow) 是计算机能表示的大于 0 的最小浮点数．下溢值除以 2 之后，在计算机中会因为舍入而得到 0．所以为了寻找下溢值，只需从 1 出发，不断除以 2，即可得到下溢值．

3．上溢值(overflow) 是计算机能表示的最大浮点数．上溢值乘以 2 之后，在计算机中会保存为 +Inf（有些可能保存为NaN）．所以为了寻找上溢值，只需从 1 出发，不断乘以 2，即可得到上溢值．

减分点：

  • 不理解概念，判断出错．

问题2

根据精度丢失定理知：当$1-\dfrac{\sin x}{x}\le \dfrac{1}{2}$时，位的丢失可以限定在1位．

所以对于不同的$x$，设法满足上述不等式．

①当$\vert x\vert \ge 2$时，有 $1-\dfrac{\sin x}{x}\le\dfrac{1}{2}$，直接调用库函数．

②当$\vert x\vert < 2$ 时，计算 $x-\sin x$ 的Taylor展开

\[\dfrac{x^3}{6}-\dfrac{x^5}{120} + \dfrac{x^7}{5040}-\cdots.\]

然后分析对于不同的 $x$，截断至多少项来进行近似计算．当 $\vert x\vert < 2$ 时，被截断的项满足

\[\begin{aligned} \left\vert\sum\limits_{n=k}^{\infty}(-1)^k\dfrac{x^{2n+1}}{(2n+1)!}\right\vert &\le \vert x \vert^{2k+1}\sum\limits_{n=k}^{\infty}\dfrac{1}{(2n+1)!} \\ &\le \dfrac{2\vert x\vert^{2k+1}}{(2k+1)!} \\ &\le \dfrac{2^{2k+2}}{(2k+1)!}. \end{aligned}\]

（也可以用带余项的Taylor展开式来得到与上式类似的不等式）

所以可以先取定 $k$ 使得 $\dfrac{2^{2k+2}}{(2k+1)!}$ 小于机器精度，这样，计算 $y=x-\sin x$ 与计算

\[p_k(x)=\sum\limits_{n=1}^{k-1}(-1)^k\dfrac{x^{2n+1}}{(2n+1)!}\]

之间的截断误差可忽略不计．

计算多项式 $p_k(x)$ 可采用秦九韶算法(Horner算法)．利用这个算法，会出现舍入误差．

用Horner算法计算 $p_k(x)$ 的步骤如下：

  • Step 1：计算 $y\leftarrow 1$；
  • Step 2：对 $n=k-1,k-2,\cdots,2$，
    计算 $y\leftarrow 1-\dfrac{x^2}{(2n)(2n+1)}\cdot y$
  • Step 3：计算 $y\leftarrow \dfrac{x^3}{6}\cdot y$

比如当 $k=4$ 时，得到的就是

\[p_4(x)=\dfrac{x^3}{6}-\dfrac{x^5}{120}+\dfrac{x^7}{5040} = \dfrac{x^3}{6}\left[1-\dfrac{x^2}{20}\left(1-\dfrac{x^2}{42}\right)\right]\]

计算 $p_k(x)$ 时，用“精度丢失定理”来判断丢失的精度不超过1位．这是因为，Step 2 中每一步浮点数运算都满足

\[1-\dfrac{x^2}{(2n)(2n+1)}\cdot y > \dfrac{1}{2},\]

最后一步也满足，所以丢失的精度不超过1位（具体写起来其实是比较繁的）．

加分点： 对舍入误差进行了理论分析的同学可以获得加分，没分析也不扣分（毕竟本题只是要求“编程验证”，没要求必须做理论分析）．

减分点：

  • 算法思路有误．
  • 没说清楚 $\vert x\vert \le 2$ 该怎么进行近似计算；
  • 绝大部分同学都没有做数值试验验证丢失的精度不超过1位（具体验证方法，请看本页面开头的“答疑”）；
  • 计算泰勒展开时，没有用 Horner 算法，直接去计算 $\dfrac{x^{2n+1}}{(2n+1)!}$．在计算机中计算阶乘是很容易崩掉的．
  • 放缩错误：当 $0 < x < 2$ 时，写 $1-\dfrac{\sin x}{x} \le \dfrac{1}{3!}x^2-\dfrac{1}{5!}x^4 \le \dfrac{1}{2}$；
  • 只是截断到 $\dfrac{x^3}{3!}-\dfrac{x^5}{5!}$ 就断言得到的结果与 $x-\sin x$ 的真实值接近．然而如果这样截断，从截断误差的角度来看根本就不接近！

没必要写的点：

  • 有同学分析了不等式 $1-\dfrac{\sin x}{x} \ge \dfrac{1}{2}$ 成立的 $x$ 的取值范围是 $\vert x\vert \ge 1.89549427$，这一步没什么必要．直接用 $x\ge 2$ 就行了．

问题3

理论上，如果初始 $\tilde{y}_0$ 有误差 $\varepsilon_0$，其满足

\[\varepsilon_0=\vert y_0-\tilde{y}_0\vert,\]

令 $\varepsilon_n = \vert y_n-\tilde{y}_n\vert$，则这部分的误差会被放大为

\[\varepsilon_n = n\varepsilon_{n-1} = \cdots = n!\varepsilon_0.\]

所以可以发现当迭代几步之后，误差就会很大．当 $n=14$ 时甚至会得到负数．

减分点： 只描述了数值结果，没有对结果的分析．

问题4

要验证算法是否数值稳定，需要对初值做一个微小的扰动（比如把 $x_0$ 改成 $x_0\pm 10^{-6}$，观察后续输出的变化情况），如果两次输出相差较大，说明这个算法不稳定；如果相差较小，说明算法稳定．

理论上，设计算机中迭代序列为 $\lbrace \tilde{x}_n\rbrace$，初始 $\tilde{x}_0$ 有误差 $\varepsilon_0$，初始 $\tilde{x}_1$ 有误差 $\varepsilon_1$，

令 $\varepsilon_n = \tilde{x}_n-x_n$，则初始部分的误差会变成

\[\varepsilon_n = \dfrac{3\varepsilon_1-\varepsilon_0}{11}\times 4^n + \dfrac{-3\varepsilon_1+12\varepsilon_0}{11}\times(\dfrac{1}{3})^n. \qquad (*)\]

当 $n\to\infty$ 时，上式等号右边第二项趋于0，初始部分的误差被缩小了．但如果 $3\varepsilon_1-\varepsilon_0\ne 0$，那么右边第一项就趋向于无穷，初始部分的误差被放大了．

从相对误差的角度去考虑： 可算出采用第一组初值的真解是$x_n=(\dfrac{1}{3})^n$，采用第二组初值的真解是$x_n=4^n$．根据$(*)$，对于第一组初值，相对误差 $r_n$ 满足

\[r_n = \dfrac{3\varepsilon_1-\varepsilon_0}{11}\times 12^n + \dfrac{-3\varepsilon_1+12\varepsilon_0}{11} ,\]

所以相对误差会迅速增长．对于第二组初值，相对误差 $r_n$ 满足

\[r_n = \dfrac{3\varepsilon_1-\varepsilon_0}{11} + \dfrac{-3\varepsilon_1+12\varepsilon_0}{11} (\dfrac{1}{12})^n ,\]

所以相对误差大约保持在 $\dfrac{3\varepsilon_1-\varepsilon_0}{11}$ ，是个非常小的值．

计算顺序的影响： 在计算 $x_{n}=\dfrac{13}{3}x_{n-1}-\dfrac{4}{3}x_{n-2}$ 时，有下面两种实现方式：

  • 方式一：按 $x_{n}=\dfrac{13x_{n-1}-4x_{n-2}}{3}$ 计算；
  • 方式二：按 $x_{n}=\dfrac{13}{3}x_{n-1}-\dfrac{4}{3}x_{n-2}$ 计算．

其 Python 代码参见附录1，采用 $x_0=1$，$x_1=4$ 的输出结果（部分）如下（完整输出结果参见附录1）：

1.0000000000000000                  	 1.0000000000000000                 
4.0000000000000000                  	 4.0000000000000000                 
16.0000000000000000                 	 15.9999999999999982                
64.0000000000000000                 	 63.9999999999999787                
256.0000000000000000                	 255.9999999999998863               
……………………………………………………………………………………………………………………（中间省略）
4503599627370496.0000000000000000   	 4503599627370494.0000000000000000  
18014398509481984.0000000000000000  	 18014398509481976.0000000000000000 
72057594037927936.0000000000000000  	 72057594037927904.0000000000000000 
288230376151711744.0000000000000000 	 288230376151711616.0000000000000000
1152921504606846976.0000000000000000 	 1152921504606846464.0000000000000000
4611686018427387904.0000000000000000 	 4611686018427385856.0000000000000000

从舍入误差的角度考虑： 在本题中，第一组初值满足 $3\varepsilon_1-\varepsilon_0\ne 0$．

第二组初值满足 $3\varepsilon_1-\varepsilon_0 = 0$．如果采用方式一，后续得到的数都是整数，所以没有舍入误差；如果采用方式二，后续得到的数会出现精度丢失，这是因为计算机中无法精确存储 $\dfrac{13}{3}$ 与 $\dfrac{4}{3}$，计算过程会带来误差．但是无论用方式一还是方式二，做微小扰动都可以让相对误差一直保持在一个很小的范围，所以我们可以认为这个算法在第二组初值下是稳定的．

因此，这个算法的稳定性与所给初值有关．

加分点： 探究了对于初值 $x_0=1$，$x_1=4$，采用方式一和方式二，说明计算精度会受到运算顺序的影响．

减分点：

  • 只描述了数值结果，没有对结果的分析．
  • 没有通过微扰初值来验证稳定性（或者没用理论说明初值对稳定性的影响），只是在验证误差很小．
  • 用第一组初值求出了$\lbrace x_n\rbrace$的通项公式 $x_n=(\dfrac{1}{3})^n$，然后就说数列递推式中的减法出现了“相减相消”现象．但是$\lbrace x_n\rbrace$并不是一直都很小的，迭代到 $n=15$ 时 $x_n$ 是一个比较大的数，就不再是相减相消了．也就是相当于没有揭示算法不稳定的根本原因．合理的解释是用舍入误差来做分析．
  • 用条件数来分析，但是该同学先导出 $\lbrace x_n\rbrace$ 的通项公式
    $x_n(t)=\dfrac{12-3t}{11}\times(\dfrac{1}{3})^n + \dfrac{3t-1}{11}\times 4^n$，
    其中 $x_0=1$，$x_1=t$，然后分析了 $x_n(t)$ 的条件数．这样的做法看似合理，然而，本题并不是直接计算 $x_n(t)$，而是计算一个数列递推式．应当针对数列递推式进行分析．

附录1：问题4的Python代码及完整输出结果

import numpy as np

## Case 1
#x0 = 1
#x1 = 1.0 / 3

## Case 2
x0 = 1
x1 = 4

y0 = x0
y1 = x1

print("{:<35.16f} \t {:<35.16f}".format(x0,y0))
print("{:<35.16f} \t {:<35.16f}".format(x1,y1))

N = 30
for i in range(N):
    x2 = (13.0 * x1 - 4.0 * x0) / 3.0
    y2 = (13.0 / 3.0) * y1 - (4.0 / 3.0) * y0
    # print("{:.16f} \t {:.16f}".format(x2,y2))
    print("{:<35.16f} \t {:<35.16f}".format(x2,y2))
    x0 = x1
    x1 = x2
    y0 = y1
    y1 = y2

完整输出结果：

1.0000000000000000                  	 1.0000000000000000                 
4.0000000000000000                  	 4.0000000000000000                 
16.0000000000000000                 	 15.9999999999999982                
64.0000000000000000                 	 63.9999999999999787                
256.0000000000000000                	 255.9999999999998863               
1024.0000000000000000               	 1023.9999999999995453              
4096.0000000000000000               	 4095.9999999999981810              
16384.0000000000000000              	 16383.9999999999927240             
65536.0000000000000000              	 65535.9999999999708962             
262144.0000000000000000             	 262143.9999999998835847            
1048576.0000000000000000            	 1048575.9999999995343387           
4194304.0000000000000000            	 4194303.9999999981373549           
16777216.0000000000000000           	 16777215.9999999925494194          
67108864.0000000000000000           	 67108863.9999999701976776          
268435456.0000000000000000          	 268435455.9999998807907104         
1073741824.0000000000000000         	 1073741823.9999995231628418        
4294967296.0000000000000000         	 4294967295.9999980926513672        
17179869184.0000000000000000        	 17179869183.9999923706054688       
68719476736.0000000000000000        	 68719476735.9999694824218750       
274877906944.0000000000000000       	 274877906943.9998779296875000      
1099511627776.0000000000000000      	 1099511627775.9995117187500000     
4398046511104.0000000000000000      	 4398046511103.9980468750000000     
17592186044416.0000000000000000     	 17592186044415.9921875000000000    
70368744177664.0000000000000000     	 70368744177663.9687500000000000    
281474976710656.0000000000000000    	 281474976710655.8750000000000000   
1125899906842624.0000000000000000   	 1125899906842623.5000000000000000  
4503599627370496.0000000000000000   	 4503599627370494.0000000000000000  
18014398509481984.0000000000000000  	 18014398509481976.0000000000000000 
72057594037927936.0000000000000000  	 72057594037927904.0000000000000000 
288230376151711744.0000000000000000 	 288230376151711616.0000000000000000
1152921504606846976.0000000000000000 	 1152921504606846464.0000000000000000
4611686018427387904.0000000000000000 	 4611686018427385856.0000000000000000