跳转至

Chapter 01 : Mathematical Preliminaries

Roundoff Errors and Computer Arithmetic

  • 截断误差（Truncation Error）：使用截断的（或者说有限的）求和来近似无穷级数的和产生的误差
  • 舍入误差（Roundoff Error）：当计算机执行实数计算时产生的误差。这是因为计算机中的算术运算涉及到的数字只有有限位数

理论上截断误差和舍入误差有一定的重合（比如 \(\frac{1}{3}=0.3333...\approx 0.3333\)），那么如何区分呢？截断误差更注重级数的“和”而舍入误差更注重数字的运算

Normalized Decimal Floating-point Form of a Real Number

对于一个 \(k\) 位十进制机器数（\(k-\) digit Decimal Machine Number）都可以被表示为 \(\pm 0.d_1d_2...d_k\times 10^n\)，其中 \(1\leq d_1\leq 9,0\leq d_i\leq 9(i=2,...,k)\)

那么对于一个实数 \(y=0.d_1d_2...d_kd_{k+1}d_{k+2}...\times 10^n\) 来说，截断法和舍入法可以表示为：

\[ fl(y)=\begin{cases} 0.d_1d_2...d_k\times 10^n(\text{Chopping})\\ \text{Chop}(y+5\times 10^{n-(k+1)})=0.\delta_1\delta_2...\delta_k\times 10^n \end{cases} \]

Absolute Error and Relative Error

  • 如果 \(p^*\) 是 \(p\) 的近似值，定义 \(|p^*-p|\) 为绝对误差（Absolute Error）， \(\frac{|p^*-p|}{|p|}(p\not=0)\) 为相对误差（Relative Error）
  • 如果有 \(t\) 是满足 \(\frac{|p^*-p|}{|p|}<5\times 10^{-t}\) 的最大非负整数，那么称 \(p^*\) 是 \(p\) 精确到 \(t\) 位有效数字（Significant Digits/Figures）的估计值

近似产生的误差

  • Chopping : \(|\frac{y-fl(y)}{y}|=|\frac{0.d_1d_2...d_kd_{k+1}...\times 10^n-0.d_1d_2...d_k\times 10^n}{0.d_1d_2...d_kd_{k+1}...\times 10^n}|=|\frac{0.d_{k+1}d_{k+2}...}{0.d_1d_2...}|\times 10^{-k}\leq\frac{1}{0.1}\times 10^{-k}=10^{-k+1}\)
  • Rounding : \(|\frac{y-fl(y)}{y}|\leq\frac{0.5\times 10^{n-k}}{0.d_1d_2...d_k\times 10^{n}}=\frac{0.5}{0.1}\times 10^{-k}=0.5\times 10^{-k+1}\)

误差的产生

  • 两个近乎相等的数相减会导致有效数字的相消，会有较大的相对误差

例如：\(a_1=0.12345,a_2=0.12346\)（这里的等号并非真实“等号”），两者本来都有 \(5\) 位有效数字，但是相减后只剩下 \(1\) 位有效数字

我们记 \(a_1=0.12345+\epsilon_1,a_2=0.12346+\epsilon_2(\epsilon_1,\epsilon_2\in [-5\times 10^{-6},5\times 10^{-6}])\)（这里的等号才是真正的等号），那么有 \(e=\epsilon_2-\epsilon_1\in [-10^{-5},10^{-5}]\)

这时候 \(a_2-a_1\) 的近似值为 \(0.00001\)，真实值为 \(0.00001+e\)，相对误差 \(e_{\text{rel}}=\frac{e}{0.00001+e}\leq 0.5\)

  • 如果有限位的表示或计算产生了误差，则除以一个小数（或者乘以一个大数）会放大绝对误差

例如，对于一个实数 \(\frac{p}{q}\)，如果 \(q\) 有一定的扰动 \(\epsilon_q\)，那么绝对误差 \(e=\frac{p}{q+\epsilon_q}-\frac{p}{q}\)，这可以近似视为 \((\frac{p}{q})'=-\frac{p}{q^2}\)，那么可以看到如果 \(q\) 比较小的话那么绝对误差会较大

误差计算举例

计算 \(f(x)=x^3-6.1x^2+3.2x+1.5\) 在 \(x=4.71\) 的值，精确到三位有效数字

\(x\) \(x^2\) \(x^3\) \(6.1x^2\) \(3.2x\)
Exact \(4.71\) \(22.1841\) \(104.487111\) \(135.32301\) \(15.072\)
Chopping \(4.71\) \(22.1\) \(104.\) \(6.1*22.1=134.81\approx 134.\) \(15.0\)
Rounding \(4.71\) \(22.2\) \(22.2*4.71=104.562\approx 105.\) \(6.1*22.2=135.42\approx 135.\) \(15.1\)
  • Exact Value : \(f(4.71)=104.487111−135.32301+15.072+1.5=−14.263899\)
  • Chopping : \(f(4.71)=104−134+15.0+1.5=−13.5\)
    • Relative Error : \(\frac{|-14.263899+13.5|}{|-14.263899|}\approx 5.35\%\)
  • Rounding : \(f(4.71)=105−135+15.1+1.5=−13.4\)
    • Relative Error : \(\frac{|−14.263899+13.4|}{|−14.263899|}\approx 6.06\%\)

可以看到，有时候舍入误差比截断误差更大，但是这样的算法误差还是太大了。

乘法较多会导致式子的误差较大，此时可以使用秦九韶算法（其实就是不断提取 \(x\)）来减少乘法的次数

  • 秦九韶算法：将多项式 \(f(x)=a_nx^n+a_{n-1}x^{n-1}+...+a_1x+a_0\) 写成形式 \(f(x)=((...(a_nx+a_{n-1})x+a_{n-2})x+...)x+a_1)x+a_0\)，然后从最内层开始计算

将例子改写成秦九韶算法：\(f(x)=((x-6.1)x+3.2)x+1.5\)

  • Chopping : \(f(4.71)=((4.71)-6.1)*4.71+3.2)*4.71+1.5=(-1.39*4.71+3.2)*4.71+1.5=(-6.54+3.2)*4.71+1.5=-3.34*4.71+1.5=-15.7+1.5=-14.2\)
    • Relative Error : \(\frac{|−14.263899+14.2|}{|−14.263899|}\approx 0.44\%\)
  • Rounding : \(f(4.71)=((4.71)-6.1)*4.71+3.2)*4.71+1.5=(-1.39*4.71+3.2)*4.71+1.5=(-6.55+3.2)*4.71+1.5=-3.35*4.71+1.5=-15.8+1.5=-14.3\)
    • Relative Error : \(\frac{|−14.263899+14.3|}{|−14.263899|}\approx 0.25\%\)

可以看到，此时误差明显减小了。

Algorithms and Convergence

Stability

  • 一个算法，如果初始数据的小变化会导致最终结果的小变化，则称为稳定（Stable）；否则称为不稳定（Unstable）
  • 一个算法，如果只有在某些初始数据的选择下才是稳定的，则称为条件稳定（Conditionally Stable）

The Growth of Errors

假设 \(E_0>0\) 是初始误差（Initial Error），\(E_n\) 是第 \(n\) 步的误差

  • 如果 \(E_n\approx CnE_0\)，其中 \(C\) 为与 \(n\) 独立的常数，则称误差的增长为线性增长（Linear Growth）
    • 线性增长的误差通常是无法避免的，当 \(C\) 和 \(E_0\) 都很小的时候，结果通常是可以接受的。
  • 如果 \(E_n\approx C^nE_0\)，其中 \(C>1\)，则称误差的增长为指数增长（Exponential Growth）
    • 指数增长的误差应该避免，因为即使 \(E_0\) 很小，\(C^n\) 也会变得很大。这会导致不可接受的不准确性。

Example

分析 \(I_n=\frac{1}{e}\int_0^1x^ne^xdx,n=0,1,2,...\)

首先我们可以考虑一下递推，我们有：

\[ I_n=\int_0^1x^nde^x=x^ne^x|_0^1-\int_0^1e^xdx^n=e-n\int_0^1x^{n-1}e^xdx=e-nI_{n-1} \]

此时 \(I_0=\frac{1}{e}\int_0^1e^xdx=1-\frac{1}{e}\approx 0.63212056=I_0^*\Rightarrow|E_0|=|I_0-I_0^*|<0.5\times 10^{-8}\)

而 \(I_1^*=1-1·I_0^*=0.36787944,...,I_{10}^*=1-10·I_9^*=0.08812800,I_{11}^*=1-11·I_{10}^*=0.03059200,I_{12}^*=1-12·I_{11}^*=0.63289600,I_{13}^*=1-13·I_{12}^*=-7.2276480,I_{14}^*=1-14·I_{13}^*=94.959424,I_{15}^*=1-15·I_{14}^*=-1423.3914\)

但是，我们却有 \(\frac{1}{e}\int_0^1x^n·e^0dx<I_n<\frac{1}{e}\int_0^1x^n·e^1dx\Rightarrow\frac{1}{e(n+1)}<I_n<\frac{1}{n+1}\)，这与我们算出的结果完全矛盾，这说明这样的算法是错误的。

这是因为 \(|E_n|=|I_n-I_n^*|=|(1-nI_{n-1})-(1-nI_{n-1}^*)|=n|E_{n-1}|=...=n!|E_0|\)，很显然，这样的误差增长是不稳定的！

由 \(I_n=1-nI_{n-1}\Rightarrow I_{n-1}=\frac{1}{n}(1-I_n)\) 以及 \(\frac{1}{e(N+1)}<I_N<\frac{1}{N+1}\)，我们倒着计算，取 \(I_N^*=\frac{1}{2}[\frac{1}{e(N+1)}+\frac{1}{N+1}]\approx I_n\)，当 \(N\rightarrow +\infty\) 时 \(|E_N|=|I_N-I_N^*|\rightarrow 0\)

取 \(I_{15}^*=\frac{1}{2}[\frac{1}{e·16}+\frac{1}{16}]\approx 0.042746233\Rightarrow I_{14}^*=\frac{1}{15}(1-I_{15}^*)\approx 0.063816918,I_{13}^*=\frac{1}{14}(1-I_{14}^*)\approx 0.066870220,I_{12}^*=\frac{1}{13}(1-I_{13}^*)\approx 0.071779214,I_{11}^*=\frac{1}{12}(1-I_{12}^*)\approx 0.077351732,I_{10}^*=\frac{1}{11}(1-I_{11}^*)\approx 0.083877115,...,I_{1}^*=\frac{1}{2}(1-I_{2}^*)\approx 0.36787944,I_{0}^*=\frac{1}{1}(1-I_{1}^*)\approx 0.63212056\)

此时我们有 \(|E_{N-1}|=|\frac{1}{N}(1-I_N)-\frac{1}{N}(1-I_N^*)|=\frac{1}{N}|E_N|\Rightarrow |E_n|=\frac{1}{N(N-1)...(n+1)}|E_N|\)，这样的误差是稳定的！

评论