一、误差的基本概念
由数学方法解决实际问题时，通常按照以下过程：
$$ 实际问题\xrightarrow{抽象、简化}数学模型\xrightarrow{数值计算}问题近似解 $$
引起误差的原因有很多：
- 模型误差：实际问题的解与数学模型解的之差。
- 观测误差：数学问题的一些参量的值往往由观测得到，但是观测不可能绝对准确，由此产生的误差称为“观测误差”。
- 截断误差：一般数学问题难以求出精确解，需要简化为较易求解的问题，以简化问题的解作为原问题解的近似，例如泰勒展开后省略后面的无穷多项。
- 舍入误差：计算过程中，受到机器字长的限制，无穷小数和位数很多的数必须舍入成一定的位数，这样的误差为“舍入误差”。
二、绝对误差和相对误差
设$x^*$为准确值$x$的一个近似值，称：
$$ e(x^*)=x-x^* $$
为近似值$x^*$的绝对误差，简称误差。
根据测量工具或计算情况，可以估计出$e(x^*)$的取值范围，即估出误差绝对值的一个上界：
$$ |e(x^*)|=|x-x^*|\le\varepsilon $$
称$\varepsilon$为近似值$x^*$的绝对误差限，绝对误差限越小越好。根据绝对误差限的定义，可以将准确值写作：
$$ x=x^*\pm\varepsilon $$
设$x^*$为准确值$x$的近似值，称绝对误差与准确值之比为近似值的相对误差，记作$e_r(x^*)$：
$$ e_r(x^*)=\frac{e(x^*)}{x}=\frac{x-x^*}{x} $$
由于计算过程中准确值$x$难以得知。所以一般取相对误差为：
$$ e_r(x^*)=\frac{e(x^*)}{x^*} $$
证明相对误差可以如上取得：
当$|e(x^*)|$很小时，可以得到：
$$ \begin{aligned} \frac{e(x^*)}{x^*}-\frac{e(x^*)}{x}&=\frac{e(x^*)(x-x^*)}{xx^*}\\&=\frac{e^2(x^*)}{xx*}\\&=e_{r}(x^*)\cdot\tilde{e}_r(x^*) \end{aligned} $$
可知两者的差是$e_r(x^*)$的高阶无穷小，可以忽略不计。
同样，相对误差只能估计其上限，如果存在正数$\varepsilon_r$使得：
$$ |e_r(x^*)|=|\frac{x-x^*}{x^*}|\le\varepsilon_r $$
称其为$x^*$的相对误差限，显然，误差限与近似值绝对值之比$\frac{\varepsilon}{|x^*|}$为$x^*$的一个相对误差限。
三、误差的传播
1、基本运算中的误差估计
由微分学，当自变量的改变量（误差）很小时，函数的微分作为函数的改变量的主要线性部分可以近似函数的改变量，故可以利用微分运算公式导出误差运算公式。设数值计算中求得的解与参量$x_1,x_2,\cdots,x_n$有关，记为：
$$ y=f(x_1,x_2,\cdots,x_n) $$
设不同参量的近似值为$x_1^*,x_2^*,\cdots,x_n^*$，相应的解也会产生误差：
$$ y^*=f(x_1^*,x_2^*,\cdots,x_n^*) $$
假设$f$在点$(x_1^*,x_2^*,\cdots,x_n^*)$处可微，则当参量误差很小时，根据：
$$ f(x_1,x_2,\cdots,x_n)=\nabla f(x_1^*,x_2^*,\cdots,x_n^*)\cdot(x_1-x_1^*,x_2-x_2^*,\cdots,x_n-x_n^*)+f(x_1^*,x_2^*,\cdots,x_n^*) $$
$$ \nabla f(x_1^*,x_2^*,\cdots,x_n^*)\cdot(x_1-x_1^*,x_2-x_2^*,\cdots,x_n-x_n^*)\approx\mathrm{d}f(x_1^*,x_2^*,\cdots,x_n^*) $$
解的绝对误差为：
$$ \begin{aligned} e(y^*)&=y-y^*\\&=f(x_1,x_2,\cdots,x_n)-f(x_1^*,x_2^*,\cdots,x_n^*)\\&\approx\mathrm{d}f(x_1^*,x_2^*,\cdots,x_n^*) \end{aligned} $$
其相对误差为：
$$ \begin{aligned} e_r(y^*)&=\frac{e(y^*)}{y^*}\\&\approx\frac{\mathrm{d}f(x_1^*,x_2^*,\cdots,x_n^*)}{f(x_1^*,x_2^*,\cdots,x_n^*)}\\&=\mathrm{d}(\ln f(x_1^*,x_2^*,\cdots,x_n^*)) \end{aligned} $$
根据上面两式，可以得到和差积商的误差公式：
$$ \begin{cases} e(x_1\pm x_2)=e(x_1)\pm e(x_2)\\ e(x_1x_2)\approx x_2e(x_1)+x_1e(x_2)\\ e\left(\dfrac{x_1}{x_2}\right)\approx \dfrac{1}{x_2}e(x_1)-\dfrac{x_1}{x_2^2}e(x_2) \end{cases} $$
$$ \begin{cases} e_r(x_1\pm x_2)=\dfrac{x_1}{x_1\pm x_2}e_r(x_1)\pm \dfrac{x_2}{x_1\pm x_2}e_r(x_2)\\ e_r(x_1x_2)\approx e_r(x_1)+e_r(x_2)\\ e_r\left(\dfrac{x_1}{x_2}\right)\approx e_r(x_1)-e_r(x_2) \end{cases} $$
这样，再根据三角不等式，可以得到：
$$ |e(x_1\pm x_2)|=|e(x_1)\pm e(x_2)|\le |e(x_1)| + |e(x_2)| $$
$$ |e_r(x_1x_2)|\approx |e_r(x_1)+e_r(x_2)|\le |e_r(x_1)|+|e_r(x_2)| $$
$$ |e_r\left(\dfrac{x_1}{x_2}\right)|\approx |e_r(x_1)-e_r(x_2)|\le|e_r(x_1)|+|e_r(x_2)| $$
2、算法的数值稳定性
计算一个数学问题，求积分制：
$$ I_n=\int_0^1\frac{x^n}{x+5}\mathrm{d}x\quad(n=1,2,\cdots) $$
注意到：
$$ I_n+5I_{n-1}=\int_0^1(\frac{x^n}{x+5}+\frac{5x^{x-1}}{x+5})\mathrm{d}x=\int_0^1x^{n-1}\mathrm{d}x=\frac{1}{n} $$
根据被积分式，知随着$n$增大，$x^n$递减，$I_n$递减且大于$0$，可以得到：
$$ I_{n+1}+5I_n=\frac{1}{n+1}\Rightarrow I_n<\frac{1}{5(n+1)} $$
$$ I_{n+1}+5I_{n}=\frac{1}{n+1}\Rightarrow 6I_n >\frac{1}{n+1}\Rightarrow I_n>\frac{1}{6(n+1)} $$
即：
$$ \frac{1}{6(n+1)}<I_n<\frac{1}{5(n+1)} $$
可以设计两种算法：
- 取$I_0=\int_0^1 \frac{1}{x+5}=\ln 1.2$按照递推公式以此计算出$I_1,I_2,\cdots$的近似值。
- 取$I_n^*\approx \frac12[\frac{1}{6(n+1)}+\frac{1}{5(n+1)}]$，按照递推公式一次计算出$I_{n-1},I_{n-2},\cdots,I_0$的近似值。
使用方法一计算$I_0\sim I_{24}$，下面给出其MATLAB代码：（代码中的I(n)代表$I_{n-1}$）
format long
I = zeros(25,1);
I(1) = log(1.2);
for n = 1:24
I(n+1) = 1/n - 5 * I(n);
end
I(1:25)然后使用方法二计算，下面给出其MATLAB代码：
format long
I = zeros(25,1);
I(25) = (1 / (6 * 25) + 1 / (5 * 25)) / 2;
for n = 24:-1:1
I(n) = (1 / n - I(n+1)) / 5;
end
I(1:25)分别比较它们的输出结果：
| | 方法一 | 方法二 |
| :-: | :----------------: | :---------------: |
| 0 | 0.182321556793955 | 0.182321556793955 |
| 1 | 0.088392216030227 | 0.088392216030227 |
| 2 | 0.058038919848865 | 0.058038919848866 |
| 3 | 0.043138734089010 | 0.043138734089005 |
| 4 | 0.034306329554950 | 0.034306329554975 |
| 5 | 0.028468352225249 | 0.028468352225126 |
| 6 | 0.024324905540419 | 0.024324905541035 |
| 7 | 0.021232615155046 | 0.021232615151969 |
| 8 | 0.018836924224770 | 0.018836924240154 |
| 9 | 0.016926489987263 | 0.016926489910342 |
| 10 | 0.015367550063686 | 0.015367550448288 |
| 11 | 0.014071340590662 | 0.014071338667650 |
| 12 | 0.012976630380023 | 0.012976639995085 |
| 13 | 0.012039925022960 | 0.012039876947652 |
| 14 | 0.011228946313773 | 0.011229186690311 |
| 15 | 0.010521935097804 | 0.010520733215111 |
| 16 | 0.009890324510982 | 0.009896333924444 |
| 17 | 0.009371906856857 | 0.009341859789542 |
| 18 | 0.008696021271271 | 0.008846256607844 |
| 19 | 0.009151472591016 | 0.008400295908150 |
| 20 | 0.004242637044922 | 0.007998520459251 |
| 21 | 0.026405862394440 | 0.007626445322793 |
| 22 | -0.086574766517654 | 0.007322318840580 |
| 23 | 0.476352093457833 | 0.006866666666667 |
| 24 | -2.3400938006225 | 0.007333333333333 |
两种方法在$n$值较小的时候，结果十分相近。但是按照方法一计算的结果，在$n=22,24$时，积分值已经小于$0$，显然是错误的，这是舍入误差在计算过程中的传播所引起的后果，设$I^*$有误差$e_0$，假设在计算过程中不产生新的舍入误差，则由递推公式，可以得知第$n$项的误差为：
$$ e_n=(-5)^ne_0 $$
即每递推一次，误差的绝对值就扩大五倍，上述计算过程中采用format long有16位有效数字，故：
$$ \varepsilon_0=\frac{1}{2}\times 10^{-16} $$
那么在第$n=21$时，误差满足：
$$ \varepsilon_{21}=5^{21}>\frac{1}{2}\times 10^{-2} $$
而$I_{21}=0.026\cdots$，已经没有任何一位有效数字。而第二种方法中，尽管$I_{24}$的取值精度不高，其误差限：
$$ \varepsilon=\frac{1}{2}(\frac{1}{125}-\frac{1}{150})\approx0.00067 $$
但是每一次递推，其误差的绝对值都在减小，到$I_0$时，满足：
$$ e_0=(-\frac{1}{5})^ne_n $$
上述事实说明，对于同一数学问题，使用的算法不同，效果也不同，我们称计算过程中舍入误差不增长的算法具有数值稳定性，否则数值就是不稳定的。例如上述问题中，方法一数值不稳定，方法二数值稳定，我们希望选择数值更加稳定的方法。
四、数值计算中应该注意的问题
1、避免两个相近的数相减
根据相对误差公式：
$$ e_r(x_1 - x_2)=\dfrac{x_1}{x_1 - x_2}e_r(x_1) - \dfrac{x_2}{x_1 - x_2}e_r(x_2) $$
可知两数之差$u=x-y$的相对误差为：
$$ e_r(u)=e_r(x-y)=\frac{xe_r(x)-ye_r(y)}{x-y}=\frac{e(x)-e(y)}{x-y} $$
当$x,y$非常接近时，$u$的相对误差很大，有效数字位数严重丢失，常常需要改变计算公式。下面是一些常用变换方法：
$$ 1-\cos x=2\sin^2\frac{x}{2},\frac{1-\cos x}{\sin x}=\frac{\sin x}{1+\cos x}\quad x\to 0 $$
$$ \sqrt{x+1}-\sqrt{x}=\frac{1}{\sqrt{x+1}+\sqrt{x}},\frac 1 x-\frac 1 {x+1}=\frac{1}{x(x+1)}\quad x较大 $$
2、避免大数吃小数的情况
计算机在进行运算过程中，首先要把参加运算的数字进行对阶，即把两数都写成绝对值小于$1$而阶码相同的数，例如$a=1000001$要改写为：
$$ a=0.1\times 10^{7}+0.0000001\times 10^{7} $$
如果计算机只能表示四位小数，那么计算出来的只有$a=0.1\times 10^7$，大数部分把小数部分吃掉了。所以在连加过程中，先把较小的数字相加，然后再加大数，这样小数累加后的进位不会先被大数吃掉。
3、避免除数的绝对值远小于被除数的绝对值
根据公式：
$$ e\left(\frac xy\right)=\frac{ye(x)-xe(y)}{y^2} $$
当$|y|\ll|x|$时，舍入误差可能增大很多。
4、要简化运算、减少运算次数，提高效率
对于计算$\ln 2$的数值，如果直接使用$\ln(1+x)$的麦克劳林级数展开，取前$n$项的和用来计算近似值，截断误差为$\frac{1}{n+1}$，如果要求误差小于$10^{-5}$，则$n\ge10^5$，要对前十万项求和，计算量很大，并且舍入误差积累使得有效数字丢失严重，如果使用级数：
$$ \ln\frac{1+x}{1-x}=2x(1+\frac{1}{3}x^2+\frac{1}{5}x^4+\cdots+\frac{x^{2n}}{2n+1}+\cdots) $$
来计算，取$x=\frac 13$，取级数的前五项即可满足要求。显然这种方法更有效。
对于计算多项式的值：
$$ P_n(x)=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0 $$
如果直接相加，那么需要做$\frac{n(n+1)}2$次乘法和$n$次加法，但是改为使用秦九韶算法：
$$ P_n(x)=a_0+x\{a_1+x[a_2+x(a_3+x(a_4+\cdots)\cdots)]\} $$
则只需要进行$n$次加法和$n$次乘法。