道长的数值分析笔记：分析误差

(一) 估算误差

　一般来说误差有可能来自于以下集中情况：

• 模型误差 Modeling Error

• 观测误差 Measurement Error

• 截断误差 Trunction Error

• 舍入误差 Roundoff Error

(1.1) 直接近似算法的误差估计

　截断误差来自于近似算法本身，例如计算式子 \(\int_0^1e^{-x^2}dx\) ，由于式子的解析不容易写出，将其泰勒展开之后，再对展开式积分，得到下列式子：

\[\begin{align} \int_0^1e^{-x^2}dx&=\int_0^1(1-x^2+\cfrac{x^4}{2!}-\cfrac{x^6}{3!}+\cfrac{x^8}{4!}-...)dx\\ &=1-\cfrac{1}{3}+\cfrac{1}{2!}\times\cfrac{1}{5}-\cfrac{1}{3!}\times\cfrac{1}{7}+\cfrac{1}{4!}\times\cfrac{1}{9}-... \end{align} \]

　取其前四项，记作 \(S_4\) ，剩余的截取部分 (Excluded terms) 记作余项(Remainder)，记作 \(R_4\) 则其截断误差则有，

\[|R_4|=\cfrac{1}{4!}\times\cfrac{1}{9}-...<\cfrac{1}{4!}\times\cfrac{1}{9}<0.005 \]

　其中每一项的分数都要按照四舍五入的方式保留三位小数，也即是说，每个舍入项的误差不超过 \(10^{-3}\) 一半，总共有两项发生舍入，那么舍入的误差最多会有 \(0.5\times 10^{-3} \times 2 = 1\times 10^{-3}=0.001\)，

\[S_4\approx1-0.333+0.1-0.025=0.743 \]

　整体的误差是由截断误差与舍入误差两个部分构成，也即是说 \(0.005+0.001=0.006\)，算得整体误差之后，通常我们会把近似算法求出的结果记作一个区间形式：

\[\int_0^1e^{-x^2}dx=0.743±0.006 \]

(1.2) 迭代近似算法的误差估计

　如果式子是直接算出结果的，那么上述的误差估计也就足够了，然而数值分析领域的大部分算法均是迭代实现的，比如迭代法求解大型线性方程组、最小二乘法拟合散点的曲线，又如牛顿迭代法等，这些算法无法一次性算出结果，因而每一步的误差都有可能传到下一步，造成误差的累积。下面回忆一个分布积分的典例，在学微积分的分部积分公式的时候，我们见过下例子 \(I_n\)， 这个式子分部积分之后会得到一个递推式，

\[\begin{align} I_n &= \cfrac{1}{e}\int_0^1x^ne^{x}dx\\ &=\frac{1}{e}(x^ne^x|^1_0-n\int_0^1x^{n-1}e^xdx)\\ &=1-nI_{n-1} \end{align} \]

　其中 \(I_0=1-\cfrac{1}{e}\) 可以直接算出来，由于积分的区间位于 \([0,1]\)，函数\(x^n\) 与 \(e^x\) 在此范围之内均为正数，因而我们可以代入积分的上下界来估算这个式子的范围，

\[\cfrac{1}{e}\int_0^1x^ne^{0}dx < I_n < \cfrac{1}{e}\int_0^1x^ne^{1}dx \\ \]

　解出其上下界之后，不难发现 \(I_n\) 是单调递减的，

\[0 < \cfrac{1}{e(n+1)} < I_n < \cfrac{1}{(n+1)} \]

　对于 \(I_0\) 保留八位小数，从小到大去算 \(I_n\) ，会发现计算机求得的结果误差极其严重。由于不同编程语言底层实现不同，可能得到的数值不太一样，但当打印若干项计算结果之后，均能发现 \(I_n\) 出现了负数，而且数值远远偏离了理论估计的上下界。由于递推式不存在截断误差，因此所有的误差均来自于舍入误差，

#include <stdio.h>
#include <math.h>

float Ix[100];
float e = 2.718281828459;

float keep(float x, int n){
    const int base = pow(10, n);
    return round(x * base) / base; // 对于给定浮点数，保留小数点位数 n
}

float I(int n){
    if(Ix[n] != 0){
        return Ix[n];
    }else{
        Ix[n] = keep(1 - n * I(n - 1), 8);
        return Ix[n];
    }
}

int main() {
    Ix[0] = keep(1 - 1 / e, 8);
    for (int i = 0; i <= 15; i++) {
        printf("i = %02d, %+15.8f\n", i, I(i));
    }
    return 0;
}

　对于第 \(n\) 步误差 \(E_n\)，

\[\begin{align} |E_n| &= |I_n - \hat{I}_n|\\ &= |(1-nI_{n-1})-(1-n\hat{I}_{n-1}) =n|E_{n-1}| =n!|E_{0}| \end{align} \]

　因而这种算法是数值不稳定的，不适合用于计算求解，为此我们可以转化从较大项递推求解较小项，\(I_{n-1}=\cfrac{1}{n}(1-I_n)\)，如此一来误差变为，

\[\begin{align} |E_{n-1}| &= |I_{n-1} - \hat{I}_{n-1}|\\ &=\cfrac{1}{n}|E_{n}| \end{align} \]

　以此类推，对于 \(N < n\) 会有：

\[E_N = \cfrac{1}{n(n-1)...(N+1)}E_n \]

　然而递推需要知道 \(I_n\) 数值才能进行，不过呢，由于 \(n\) 较大的时候，其上下界已经非常接近了，所以我们可以直接取其中点作为 \(I_n\) 数值，然后反推每一项，通过这种算法得出的结果是数值稳定的，

#include <stdio.h>
#include <math.h>

float Ix[100];
float e = 2.718281828459;

float keep(float x, int n){
    const int base = pow(10, n);
    return round(x * base) / base; // 对于给定浮点数，保留小数点位数 n
}

float I(int n){
    if(Ix[n] != 0){
        return Ix[n];
    }else{
        Ix[n] = keep((1 - I(n + 1)) / n, 8); 
        return Ix[n]; // 第 0 项需要特判，以免出现除以零的情况
    }
}

float approx(int n){
    Ix[n] = ((1 / (e * (n + 1))) + (1 / (n + 1)) ) / 2;
    return Ix[n]; // 根据上下界估算第 n 项数值等于多少，只适用于 n 较大的情况
}

int main() {
    Ix[0] = keep(1 - 1 / e, 8);
    Ix[90] = approx(90);
    for (int i = 0; i <= 15; i++) {
        printf("i = %02d, %+15.8f\n", i, I(i));
    }
    return 0;
}

(1.3) 相对误差限

　由于我们一般无从得知精确值，所以绝对误差\(e(x)=|\hat{x}-x|\)很难算出来，再一个就是不同事物的数量级是不一样的，所以我们需要无量纲化处理，得到一个相对误差 \(e_r(x)=e(x)/|x|\)，其中分母常用估计值 \(\hat{x}\) 代表真实值 \(x\)，在实际的工程问题之中，对于除不尽的数据我们需要保留若干位数，然而保留多少位才算合适是需要计算分析的。我们通常希望把误差控制于某个范围之内，因而我们会指定一个最大可容忍的误差限度，简称误差限，记作 \(\epsilon\)，那么最大可容忍的相对误差即为相对误差限，

\[\epsilon_r = |\cfrac{\epsilon}{\hat{x}}| \]

　知道了相对误差限之后，我们能够估算出至少需要保留几位小数。
　如果一个近似值 \(\hat{x}\) 从第一个非零数字数起共有 \(n\) 位数字，且误差限不超过末尾的单个单位则称 \(\hat{x}\) 具有 \(n\) 位有效数字(significantw digits)。

(1.4) 绝对条件数

　对于\(A=f(x)\)，我们可以根据 \(e(x)\) 估算 \(e(A)\)，根据中值定理展开可知，

\[\begin{align} e(A) & =f'(\hat{x})(x-\hat{x})\\ &= f'(\hat{x})e(x)\\ \end{align} \]

　式子左右两边同时除以 \(f(\hat{x})\) ，然后两边凑项，配成相对误差的形式，则有，

\[\begin{align} \cfrac{e(A)}{f(\hat{x})} &= f'(\hat{x})(x-\hat{x})\\ e_r(A)&=\cfrac{f(\hat{x})-f(x)}{\hat{x} - x}\times\cfrac{\hat{x}}{f(\hat{x})}\times\cfrac{\hat{x}-x}{\hat{x}} e_r(A)=|\cfrac{f'(x)}{f(x)}x| e(x) \end{align} \]

　其中 \(|\cfrac{f'(x)}{f(x)}x|\) 称为(绝对)条件数，可以看出一个函数的绝对条件数越大，对于输入的 \(x\) 也就越敏感。

(二) 机器浮点数带来的误差

(2.1) 二进制表示

　计算机使用规格化浮点数表示实数，也就是说，计算机无法包含整个实数集，其存储的是实际值的一个近似值，对于一个给定的实数\(x\)，计算机内部存储了，

\[x\approx ± q \times 2^n \]

　其中 \(q\) 称为尾数，整数 \(n\) 称为阶码，这是一个有限的二进制数，满足不等式 \(\cfrac{1}{2}\leq q < 1\)，其展开形式可以写成，

\[0.d_1d_2d_3...d_k \times 2^n \]

　其中 \(d_1=1\)，这一位是固定的。二进制个数受到了 \(q\) 与 \(n\) 限制，比如给定\(0.d_1d_2d_3d_4 \times 2^n\)，\(n\in\{-3,-2,-1,0,1,2,3,4\}\)，那么尾数和阶码各有八种选择，根据乘法原理，我们可以得到六十四个浮点数。尾数的表示尾数的长度决定了浮点数能够准确表达的范围。

　例如，32bit单精度浮点数，尾数位数等于 24，其具有六位的十进制精度，因为 \(2^{-(24-1)} \approx 1.2 \times 10^{-7}\)，类似的，48bit 单精度浮点数尾数 40，\(2^{-(40-1)} \approx 1.8 \times 10^{-12}\)，具有十一位的十进制精度，64bit 双精度浮点数尾数 53，\(2^{-(53-1)} \approx 2.2 \times 10^{-16}\)，具有十六位的十进制数值精度。

(2.2) 减少截断误差

　由于机器浮点数经常会出现大数吃小数的情况，设有\(a=123456, b=2.189\)，由于float的精度只有小数点后六七位左右，所以必然得不到精确解 123458.189，其中9 是必然会截掉的，做一个加法产生的误差不少于0.009，做一千次这样的加法，误差达到9，

#include<stdio.h>
float a = 123456, b = 2.189;
float x[1000], sumv = 0;

int main() {
    printf("%.8f\n", a + b);
    sumv = x[0] = a;
    for (int i = 1; i < 1000; i++) {
        x[i] = b;
        sumv += x[i];
    }
    printf("%.8f\n", sumv); // 会发现通常上述算法计算得到的结果与真实值相差较大
    return 0;
}

　故在对多个浮点数求和的时候，我们可以使用Kahan算法(Kahan's Summation Formula) 计算，这种算法思想是用一个 carry 变量对上一轮运算出现的误差进行记录，并对本轮运算之前进行修正。

#include<stdio.h>
float a = 123456, b = 2.189;
float x[1000], sumv = 0, carry = 0; // 使用 carry 对丢失的低阶位进行运行补偿

int main() {
    printf("%.8f\n", a + b);
    sumv = x[0] = a;
    for (int i = 1; i < 1000; i++) {
        x[i] = b;
        float y = x[i] - carry; // 对上一轮运算造成的误差进行补偿(修正)
        float t = sumv + y;     // 如果sum很大，y很小，那么y低阶数字丢失了
        carry = (t - sumv) - y; // 此时做差得到的 t 不等于零，而是 y 丢失的低阶数字部分
        sumv = t;
      
    }
    printf("%.8f\n", sumv);
    return 0;
}

支持作者