误差分析基础

误差的来源

  • 计算前的误差
    • 模型误差：问题建模不准确产生的误差，例如动力学计算中忽略空气阻力、地球表面积计算将地球近似为球体等
    • 数据误差：常数、测量值、上一步计算结果等存在的误差
  • 计算中的误差
    • 截断误差：计算方法上产生的误差，例如使用3阶Taylor公式近似计算f(x)=sin⁡xf(x)=\sin x，4阶及以上的小量构成了截断误差
    • 舍入误差：计算时用于表示数据的位数有限

误差及其分类

设某数据的准确值为xx，计算中采取或计算出的近似值为x^\hat x，则绝对误差定义为：

e(x^)=x^−xe(\hat x)=\hat x-x

相对误差定义为：

er(x^)=x^−xxe_r(\hat x)=\frac{\hat x-x}{x}

实际情况下，准确值往往未知，此时精确的误差无法得知，往往使用估计的最大误差来讨论，称为误差限。绝对误差限和相对误差限分别记为ε(x^)\varepsilon(\hat x)、εr(x^)\varepsilon_r(\hat x)

定理：有效数字与相对误差

设对xx保留pp位有效数字得到x^\hat x，则x^\hat x的相对误差有界：

∣εr(x^)∣≤5d0×10−p|\varepsilon_r(\hat x)|\leq\frac{5}{d_0}\times 10^{-p}

其中d0d_0为第一位有效数字

定理：上述定理的逆命题

设xx的第一位有效数字为d0d_0，若近似值x^\hat x的相对误差满足

∣εr(x^)∣≤5d0+1×10−p|\varepsilon_r(\hat x)|\leq\frac{5}{d_0+1}\times 10^{-p}

则xx与x^\hat x至少满足以下两种情形之一：

  • 二者前pp位有效数字相同
  • 或二者保留pp位有效数字后的结果相同

上面的两种情形其实是对“两个数在pp位有效数字精度下相同”这句模糊表述的两种解释，例如对x=9.1243x=9.1243，x^1=9.1247\hat x_1=9.1247和x^2=9.1238\hat x_2=9.1238分别是上面的两种情形

数据传递误差与计算误差

例如，用x^\hat x代替xx近似计算f(x)f(x)，计算结果为f^(x^)\hat f(\hat x)，则误差分为两部分：

ε(f^(x^))=[f^(x^)−f(x^)]+[f(x^)−f(x)]\varepsilon(\hat f(\hat x))=\left[\hat f(\hat x)-f(\hat x)\right]+\left[f(\hat x)-f(x)\right]

第一部分为 实际计算的函数与精确函数ff不同带来的误差 ，即计算误差

第二部分为 使用的数据并非精确值带来的误差 ，即数据传递误差

计算误差的两部分：截断误差和舍入误差

考虑以下问题：使用差商近似一阶导数，即

f′(x)≈f(x+h)−f(x)hf'(x)\approx \frac{f(x+h)-f(x)}{h}

由于

f(x+h)=f(x)+hf′(x)+h22f′′(ξ), ξ∈(x,x+h)f(x+h)=f(x)+hf'(x)+\frac{h^2}{2}f''(\xi),\ \xi\in(x,x+h)

于是

f′(x)−f(x+h)−f(x)h=−h2f′′(ξ)f'(x)-\frac{f(x+h)-f(x)}{h}=-\frac{h}{2}f''(\xi)

设∣f′′(ξ)∣≤M|f''(\xi)|\leq M，则

∣f′(x)−f(x+h)−f(x)h∣≤Mh2\left|f'(x)-\frac{f(x+h)-f(x)}{h}\right|\leq\frac{Mh}{2}

这部分是截断误差，即展开式阶数不够带来的误差

εT=Mh2\varepsilon_T=\frac{Mh}{2}

每次计算函数ff会带来ϵ\epsilon的舍入误差，则总舍入误差

εR=2ϵh\varepsilon_R=\frac{2\epsilon}{h}

于是

ε=εT+εR=Mh2+2ϵh\varepsilon=\varepsilon_T+\varepsilon_R=\frac{Mh}{2}+\frac{2\epsilon}{h}

则可以取使误差最小的hh作为最终模拟步长

hopt=2ϵMh_{opt}=2\sqrt\frac\epsilon{M}

此时的误差为

εmin=2Mϵ\varepsilon_{min}=2\sqrt\frac M{\epsilon}

问题的敏感性

考虑输入数据的扰动（即数据传递误差）对问题结果的影响，若影响很大，则问题敏感或病态，否则称问题不敏感或良态；问题敏感性相对的概念为稳定性

问题的敏感性可以用条件数来反映：

cond=∣Δf/f∣∣Δx/x∣{\rm cond}=\frac{|\Delta f/f|}{|\Delta x/x|}

当Δx→0\Delta x\to 0时

cond=∣xf′(x)f(x)∣{\rm cond}=\left|\frac{xf'(x)}{f(x)}\right|

条件数越大，问题越病态

类似地可以定义绝对条件数：

condA=∣Δf∣∣Δx∣{\rm cond}_A=\frac{|\Delta f|}{|\Delta x|}

注：条件数反应的是问题的特性，与计算方法无关，因此不考虑截断误差f^(x^)−f(x^)\hat f(\hat x)-f(\hat x)

对于简单的运算可以估算数据传递误差限：

设对某多元函数y=f(x1,x2,⋯ ,xn)y=f(x_1,x_2,\cdots,x_n)，输入数据扰动得到的近似值为y^=f(x^1,x^2,⋯ ,x^n)\hat y=f(\hat x_1,\hat x_2,\cdots,\hat x_n)

注意：这里同样不考虑截断误差

则

y^−y=f(x^1,x^2,⋯ ,x^n)−f(x1,x2,⋯ ,xn)\hat y-y=f(\hat x_1,\hat x_2,\cdots, \hat x_n)-f(x_1,x_2,\cdots,x_n)

≈∑k=1n(x^k−xk)∂f∂xk(x^1,x^2,⋯ ,x^n)\approx \sum\limits_{k=1}^n (\hat x_k-x_k)\frac{\partial f}{\partial x_k}(\hat x_1,\hat x_2,\cdots,\hat x_n)

由于

∣∑k=1n(x^k−xk)∂f∂xk(x^1,x^2,⋯ ,x^n)∣≤∑k=1n∣(x^k−xk)∂f∂xk(x^1,x^2,⋯ ,x^n)∣\left|\sum\limits_{k=1}^n (\hat x_k-x_k)\frac{\partial f}{\partial x_k}(\hat x_1,\hat x_2,\cdots,\hat x_n)\right|\leq\sum\limits_{k=1}^n\left|(\hat x_k-x_k)\frac{\partial f}{\partial x_k}(\hat x_1,\hat x_2,\cdots,\hat x_n)\right|

≤ε(x^)∑k=1n∣∂f∂xk(x^1,x^2,⋯ ,x^n)∣\leq \varepsilon(\hat x)\sum\limits_{k=1}^{n}\left|\frac{\partial f}{\partial x_k}(\hat x_1,\hat x_2,\cdots, \hat x_n)\right|

于是

ε(y^)=ε(x^)∑k=1n∣∂f∂xk(x^1,x^2,⋯ ,x^n)∣\varepsilon(\hat y)=\varepsilon(\hat x)\sum\limits_{k=1}^{n}\left|\frac{\partial f}{\partial x_k}(\hat x_1,\hat x_2,\cdots, \hat x_n)\right|

这里偏导数的自变量取近似值，是因为在实际应用中近似值可知而精确值不可知

例1：长度为100的向量求和

算法1：按存储顺序求和

算法2：先按绝对值从小到大排序，再按此顺序求和

假设所求和的数为2.0,0.01,0.01,⋯ ,0.012.0,0.01,0.01,\cdots, 0.01，求和时的数据精度为2位十进制，则算法1结果为2.02.0，算法2结果为3.03.0

于是算法2更稳定

上面分析的思路是：计算数据精度越低时，数据传递误差越大，因此如果一种算法在较低精度时偏差较大，说明其不稳定

例2：计算ϕ=5−12\phi=\frac{\sqrt 5-1}{2}的前nn次幂

设n=20n=20，采取近似值x=0.618034x=0.618034，使用双精度浮点数计算，可以忽略舍入误差

算法1：直接计算xnx^n

算法2：利用递推式

{ϕn+1=ϕn−1−ϕnϕ0=1, ϕ1≈x\begin{cases}\phi^{n+1}=\phi^{n-1}-\phi^n\\\phi^0=1,\ \phi^1\approx x\end{cases}

设e(ϕn^)=ene(\hat{\phi^n})=e_n，则在算法2中

e0=0e_0=0

e1=x−ϕe_1=x-\phi

e2=−e1e_2=-e_1

e3=e1−e2=2e1e_3=e_1-e_2=2e_1

e4=e2−e3=−3e1e_4=e_2-e_3=-3e_1

类推可知

e20=−F20e1e_{20}=-F_{20}e_1

其中F20F_{20}为Fibonacci数列第20项，此时相对误差已大于11

对于算法1，

e20=x20−ϕ20=(ϕ+e1)20−ϕ20≈20ϕ19e1e_{20}=x^{20}-\phi^{20}=\left(\phi+e_1\right)^{20}-\phi^{20}\approx 20\phi^{19}e_1

由于ϕ<1\phi<1，显然

20ϕ19<<F2020\phi^{19}<<F_{20}

结论：较快速的算法不一定稳定，实际应用中需要折中考虑

算法稳定性的向后误差分析

现在用某种方法来近似函数求值问题y=f(x)y=f(x)，设得到的近似值为（仅考虑方法误差）

y^=f^(x)\hat y=\hat f(x)

此时找到一个与上述方法误差等效的数据传递误差，即找到x^\hat x使

y^=f(x^)\hat y=f(\hat x)

则

Δx=x^−x\Delta x=\hat x-x

称为方法的向后误差，向后误差越大，说明算法可以容忍的数据传递误差越大，即算法越稳定

浮点数

浮点数存储中，尾数的有效数字位数称为精度位数pp，注意由于第一位有效数字并不存储，实际精度位数为存储尾数位数+1

浮点数的机器精度为εmach=2−p\varepsilon_{\rm mach}=2^{-p}（这实际上是一个相对误差）

设指数的下限值为LL，上限值为UU，若不考虑非规范数，则下溢值为2L2^L，上溢值为(2−2−(p−1))2U(2-2^{-(p-1)})2^U

浮点数系统由进制、精度位数、指数下限、指数上限确定

定理（浮点数的相对误差限）

设xx的浮点数表示为fl(x){\rm fl}(x)，则

∣fl(x)−xx∣≤εmach\left|\frac{\mathrm{fl}(x)-x}{x}\right|\leq\varepsilon_{\rm mach}

定理（大数吃小数）

对两个浮点数x1,x2x_1,x_2，若

∣x2x1∣≤12εmach\left|\frac{x_2}{x_1}\right|\leq \frac12\varepsilon_{\rm mach}

则x2x_2对x1+x2x_1+x_2的浮点数运算值毫无影响

浮点数的抵消现象

两个值相近、有效数字位数较高的浮点数相减，得到的数准确的有效数字位数可能很低，相对误差被放大

上述现象称为浮点数的抵消

抵消现象可以通过调整算法来避免，例如求一元二次方程的根，但其中有一个根很接近00（∣4ac∣<<b2, b>0|4ac|<<b^2,\ b>0）

此时正常应用求根公式求解靠近00的根：

x1=−b+b2−4ac2ax_1=\frac{-b+\sqrt{b^2-4ac}}{2a}

则分子会发生抵消。可以通过分子有理化来解除这个抵消：

x1=2c−b−b2−4acx_1=\frac{2c}{-b-\sqrt{b^2-4ac}}

减小舍入误差的几点建议

  • 使用双精度浮点数
  • 对于运算次数很多的算法，可以采取以下建议：
    • 避免乘除法中间结果溢出
    • 避免加减法大数吃小数
    • 避免同号相近数相减（抵消）
    • 尽可能减少运算次数