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

(一) 估算误差

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

  • 模型误差 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;
}



支持作者

posted @ 2022-11-22 17:32  道长陈牧宇  阅读(479)  评论(0)    收藏  举报