计算方法笔记
基本概念
计算方法
数值计算方法是科学计算的核心，研究对象是那些在理论上有解，而无求解公式或计算量过大必须借助计算机实现的数学问题。
简单而言就是容忍一定的误差，通过一个改进的在计算机上迭代的计算方法获取原解的近似解。
而数值计算方法主要有两大方向：
- 数值算法的构造：针对数值问题研究可在计算机上执行且行之有效的计算公式（计算方法）
- 算法的理论分析：主要包括误差分析（数值问题的性态，数值方法的截断误差，舍入误差，稳定性，收敛性等）和复杂性分析
也就是说，我们希望找到一种合理的算法能够降低其复杂度的同时保证结果的合理性。
误差
首先明确，数值方法得到的都是近似解而非真解。
那么误差的定义就是真解与近似解之间的偏差，即 \[ e=x-x^* \] 其中 \(x^*\) 是近似解的表示符号。
误差分为四类：
- 模型误差（model error）
- 观测误差（observational error）
- 截断误差（truncation error）
- 舍入误差（round-off error）
前两者是客观存在的误差，这是原问题模型本身或测量方法所存在的误差；
后两者是计算方法所引起的误差，这是我们所需要关注的。
由简化问题所引起的误差称为截断误差或方法误差，例如我们通常会将三角函数进行泰勒展开 \[ \cos x=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+\cdots \] 随后截断一部分做近似计算 \[ \cos x\approx 1-\frac{x^2}{2!} \] 这样便会产生误差 \[ e=\frac{x^4}{4!}-\frac{x^6}{6!}+\cdots \] 通过交错级数审敛法可以得知这里的截断误差 \(e<\frac{x^4}{4!}\)。
舍入误差主要是针对计算机的存储问题，比如说标准通用 CPU 都支持的 IEEE 浮点数标准，可能会使得每次计算时都存在一定的舍入。
少量的舍入误差还好，但是一般数值方法都需要大量迭代，多次的舍入误差会导致结果出现很大的波动。
绝对误差与相对误差
假定 \(x^*\) 是准确值 \(x\) 的一个近似值。
那么其绝对误差（absolute error）为 \[ e=x-x^* \] 但是往往在实际计算中，我们无法确定准确值 \(x\) 的大小，也就是说我们无法确定绝对误差 \(e\)。
所以为了方便研究分析，实际会使用绝对误差限（absolute error bound）来断定其上限，符号为 \(\epsilon\)，定义为 \[ |e|=|x-x^*|\leq \epsilon \]
或者可以写作 \[ x^*-\epsilon\leq x\leq x^*+\epsilon \]
但是绝对误差不能满足所有的计算需要，因为近似值和真值可能十分小，那么得到的误差就会极其小，无法通过绝对误差的大小刻画近似值的准确程度，那么就需要相对误差（relative error） \[ \begin{align} e_r=\frac{e}{x}=\frac{x-x^*}{x}\\ e_r=\frac{e}{x^*}=\frac{x-x^*}{x^*} \end{align} \] 即绝对误差与真值或近似值的比值。
当 \(|e_r|\) 趋近于零时， \[ \frac{e}{x}-\frac{e}{x^*}=\frac{-e^2}{xx^*} \] 是 \(e_r\) 的高阶无穷小，可以忽略不计，所以无论哪个公式都是合理的，但一般默认是绝对误差与近似值的比值，因为这可计算。
同样的，也有相对误差限（relative error bound），符号为 \(\epsilon_r\)，定义为 \[ |e_r|=|\frac{e}{x^*}|\leq \epsilon_r \] 显然，误差限（error bound）不是唯一的，且 \(\frac{\epsilon}{|x^*|}\) 是 \(x^*\) 的一个相对误差限。
误差（error）函数写为 \[ e(x^*)=x-x^* \] 误差限（error bound）函数写为 \[ e(x^*)\leq\epsilon(x^*) \]
有效数字
从左到右，从第一个非零数字开始数，一共有几位就是几位有效数字
比如说 0.023 的有效数字是两位（23），而 0.0230 的有效数字是三位（230）。
对于有效数字是有舍入的（四舍五入），那么其误差限显然是小数点后最后一位的半个单位，一般写作 \(0.5\times 10^{-n}\)，这里的 \(n\) 不是上图的 \(n\) 个数，而是小数点后的位数。
比如说 0.023 的误差限是 0.0005，0.0230 的误差限是 0.00005。
有效数字（对于近似值而言）是有严格定义的，一定是满足四舍五入的舍入规则才能被作为有效数字，以 \(x=\sqrt{3}=1.7320508\cdots\) 为例，如果说有三个近似值 \(x_1=1.73\)，\(x_2=1.7321\)，\(x_3=1.7320\)，可以计算这三个近似值的误差限分别为 \[ \begin{align} |x-x_1|=|\sqrt{3}-1.73|=0.20508\cdots\times10^{-2}\leq0.5\times10^{-2}\\ |x-x_2|=|\sqrt{3}-1.7321|=0.508\cdots\times10^{-4}\leq0.5\times10^{-4}\\ |x-x_3|=|\sqrt{3}-1.7320|=0.508\cdots\times10^{-4}\leq0.5\times10^{-3}\\ \end{align} \] 故对于三个近似值而言，分别是
- \(x_1\) 精确到小数点后 2 位，有 3 位有效数字；
- \(x_2\) 精确到小数点后 4 位，有 5 位有效数字；
- \(x_3\) 精确到小数点后 3 位，有 4 位有效数字。
误差的传播
自变量的误差势必会传播到因变量，比如说 \(y=ax+b\)，假定自变量存在误差 \(e=x-x^*\)，那么有 \[ y-y^*=a(x-x^*) \] 即 \(e(y^*)=ae(x^*)\)，对于一个线性函数，误差会直接扩大 \(a\) 倍。
扩展到对一个连续可导函数 \(f(x)\)，有 $$ \[\begin{align} \Delta y=y-y^*=f(x)-f(x^*)\\ \end{align}\] \[ 当 $\Delta y$ 很小时，可以做一个线性近似得到 \] yf'(x)(x-x^)=f'(x)x \[ 也就是 \] e(y)f'(x^)e(x) \[ 那么其误差限也应该有 \] (y)|f'(x^)|(x) \[ 而对于相对误差，有 \] e_r(y^)=y|{y=y} \[ 当 $f$ 是多元函数时，绝对误差为 \] e(y^)=n{k=1}e(x^*) \[ 相对误差为 \] \[\begin{align} e_r(y^*)&=\frac{e(y^*)}{y^*}\\ &=\sum^n_{k=1}\frac{\partial f(x_1^*,\cdots,x_n^*)}{\partial x_k}\frac{e(x^*)}{y^*}\\ &=\sum^n_{k=1}\frac{\partial f(x_1^*,\cdots,x_n^*)}{\partial x_k}\frac{x^*}{y^*}e_r(x^*) \end{align}\] $$
四则运算的误差传播
加减法的误差传播公式 \[ \begin{cases} e(x_1\pm x_2)=e(x_1)\pm e(x_2)\\ e_r(x_1\pm x_2)=\frac{x_1}{x_1\pm x_2}e_r(x_1)\pm \frac{x_2}{x_1\pm x_2}e_r(x_2)\\ \end{cases} \] 如果扩展到 \(n\) 元，显然满足叠加定理，且相对误差按权重比例分配误差影响。
乘除法的误差传播公式 \[ \begin{cases} e(x_1x_2)\approx x_2e(x_1)\pm x_1e(x_2)\\ e_r(x_1x_2)\approx e_r(x_1)+e_r(x_2)\\ \end{cases} \]
\[ \begin{cases} e(\frac{x_1}{x_2})\approx \frac{1}{x_2}e(x_1)-\frac{x_1}{x_2^2}e(x_2)\\ e_r(x_1x_2)\approx e_r(x_1)-e_r(x_2)\\ \end{cases} \]
之所以乘除法的误差传播公式是约等于，是因为误差的传播公式是做线性近似，而这里的函数是二次的，如果对于高次函数（或者说多元乘法）误差限会很大。
误差传播例题
1
已知某物体行程 \(s\) 的近似值 \(s^*=800\mathrm{m}\)，所需时间 \(t\) 的近似值 \(t^*=35\mathrm{s}\)。如果已知 \(|s-s^*|\leq 0.5\)，\(|t-t^*|\leq 0.05\)，试求平均速度 \(v\) 的绝对误差限和相对误差限。
提取题目信息，\(\epsilon(s^*)=0.5,\epsilon(t^*)=0.05,\epsilon_r(s^*)=0.5/800,\epsilon_r(t^*)=0.05/35\)
首先计算平均速度 \(v^*\) 的绝对误差 \[ \begin{align} v^*&=\frac{s^*}{t^*}\\ e(v^*)&=e(\frac{s^*}{t^*})\\ &\approx \frac{1}{t^*}e(s^*)-\frac{s^*}{t^{*2}}e(t^*) \end{align} \] 根据绝对值不等式，有 \[ \begin{align} |e(v^*)|&\leq \frac{1}{t^*}e(s^*)+\frac{s^*}{t^{*2}}e(t^*)\\ &\leq \frac{1}{t^*}\epsilon(s^*)+\frac{s^*}{t^{*2}}\epsilon(t^*)\approx 0.0469\\ \end{align} \] 即 \(\epsilon(v^*)=0.0469\)。
那么 \[ \begin{align} \epsilon_r(v^*)=\frac{\epsilon(v^*)}{v^*}=0.00205 \end{align} \]
或者也可以套用公式 \[ \begin{align} e_r(v^*)&= e_r(\frac{s^*}{t^*})\\ &\approx e_r(s^*)-e_r(t^*) \end{align} \] 那么根据绝对值不等式，有 \[ \begin{align} |e_r(v^*)|&\leq |e_r(s^*)|+|e_r(t^*)|=0.00205 \end{align} \] 即 \(\epsilon_r(v^*)=0.00205\)。
2
已知 \(a=1.2031,b=0.978\) 是四舍五入的近似值，求 \(a+b\) 和 \(a\times b\) 各有几位有效数字？
根据有效数字的定义，得知 \[ \epsilon(a)=0.5\times 10^{-4}\\ \epsilon(b)=0.5\times 10^{-3}\\ \] 那么有 \[ \begin{align} |e(a+b)|=|e(a)+e(b)|\leq \epsilon(a)+\epsilon(b)=0.55\times 10^{-3}\leq 0.5\times 10^{-2}\\ |e(a\times b)|=|be(a)+ae(b)|\leq b\epsilon(a)+a\epsilon(b)=0.65045\times 10^{-3}\leq 0.5\times 10^{-2}\\ \end{align} \] 故全都精确到小数点后 2 位，大约计算 \(1.20+0.97>2,1.20\times 0.97>1\)，故结果均为三位有效数字。
3
设计算球体积要使相对误差限为 1%，问度量半径为 \(R\) 时允许的相对误差限是多少？
首先明确球体积公式为 \[ V(R)=\frac{4}{3}\pi R^3 \] 那么绝对误差为 \[ e[V(R)]=\frac{4}{3}\pi(3R^2)e(R)=4\pi R^2e(R) \] 那么相对误差为 \[ e_r[V(R)]=\frac{e[V(R)]}{V(R)}=3\frac{e(R)}{R}=3e_r(R) \] 那么 \[ e_r(R)=\frac{1}{3}e_r[V(R)]\leq \frac{1}{3}\epsilon[V(R)]=0.3333\% \]
误差改善
为了减少误差的影响，设计算法时需要考虑以下几个方面
- 选择数值稳定的算法
- 避免两个相近的数相减
- 避免大数吃小数现象
- 避免绝对值太小的数做除数
- 简化计算步骤，减少运算次数
选择数值稳定的算法
一个算法，如果在一定的条件下，其舍入误差在整个运算过程中能够得到控制或者舍入误差的增长不影响产生可靠的结果，则称该算法是数值稳定的，否则称为数值不稳定的。
比如说计算积分式的值 \[ I_n=\int^1_0 x^ne^{x-1}\mathrm{d}x \] 利用分部积分法可以求得递推公式 \[ I_n=1-nI_{n-1} \] 其误差为 \[ e(I_n^*)=-ne(I_{n-1}^*) \] 可以尝试数值实验：
不难发现在 \(n=8,9\) 时，误差就变得难以接受了。
显然分析可得，误差传播到 \(n\) 时， \[ |e(I_n^*)|=n!|e(I_0)^*| \] 以阶乘倍扩大误差，算法不稳定。
但是反之，从高积分值求得低积分值，即递推公式为 \[
I_{n-1}=\frac{1-I_n}{n}
\] 可知误差为 \[
e(I_{n-1}^*)=-\frac{e(I_n^*)}{n}
\] 
那么误差传播到 \(0\) 时， \[ |e(I_0^*)|=\frac{|e(I_n)^*|}{n!} \] 误差反而会越来越小。
对于同一数学问题，算法设计至关重要。只有选用数值稳定性好的算法，舍入误差才能够得到控制，求得比较精确的结果。
与数值稳定性相关的概念是良态和病态问题。
对于一个数值问题，如果输入数据的细微变化(或误差)，引起输出数据(即所求解)的很大变化，则称之为病态问题；如果输入数据的细微变化，不会引起解的很大变化，则称之为良态问题。
避免两个相近的数相减
两数差的相对误差关系式为 \[ e_r(x-y)=\frac{e(x)-e(y)}{x-y} \] 如果 \(x-y\) 趋近于 0，那么相对误差将会很大，有效数字位数将严重丢失。
一般来说会用三种方式防止这种现象产生：
因式分解 \[ \begin{align} \frac{1}{x}-\frac{1}{1+x}=\frac{1}{x(x+1)} \end{align} \]
分子分母有理化 \[ \begin{align} \sqrt{x+1}-\sqrt{x}=\frac{1}{\sqrt{x+1}+\sqrt{x}}\\ \end{align} \]
三角恒等式 \[ \begin{align} 1-\cos{x}=2\sin^2\frac{x}{2}\\ \frac{1-\cos{x}}{\sin x}=\frac{\sin x}{1-\cos{x}}\\ \end{align} \]
其他 \[ \ln{x_1}-\ln{x_2}=\ln{\frac{x_1}{x_2}} \]
例如说求二次方程 \(x^2-16x+1=0\) 的根，取 \(\sqrt{63}\approx 7.94\)，那么根据求根公式有 \[ x_1=8+\sqrt{63} \] 而第二个根为 \[ x_2=8-\sqrt{63}\approx 0.06 \] 但是如果避免减法 \[ x_2=\frac{1}{8+\sqrt{63}}\approx\frac{1}{15.94}\approx 0.0627 \] 反而精度更高。
避免大数“吃”小数现象
比如说计算机实例，使用 JavaScript
1 | 0.1+0.00000000000000001 |
结果令人吃惊
避免绝对值太小的数作除数
\[ e(\frac{x}{y})\approx \frac{1}{y}e(x)-\frac{x}{y^2}e(y)\\ \]
当 \(x\ll y\) 或 \(y\to 0\) 时，\(e(\frac{x}{y})\) 可能很大。
简化计算步骤，减少运算次数
比如说快速幂算法，本身求解 \(x^{31}\) 需要 30 次乘法运算，但是使用快速幂算法仅需要 8 次乘法运算 \[ x^{31}=x^1x^2x^4x^8x^{16} \]
插值法
插值问题
插值法是数值计算的基本课题，目标是使用简单函数为离散数组建立连续模型（连续关系）。
比如说在物理实验中测得一组被测对象的关系值，比如说压力传感器的力的大小和输出电压的大小，这是一组离散点，可能数据如下
| Force (N) | voltage (mV) |
|---|---|
| 20 | 810 |
| 21 | 892 |
| 24 | 1162 |
| 25 | 1260 |
| 27 | 1468 |
| 34 | 2322 |
| 50 | 5010 |
| 52 | 5418 |
| 55 | 6060 |
| 56 | 6282 |
这时使用简单的线性关系是无法完美描述该传感器的传感特性的，使用最小二乘法做拟合得到的结果如下
其相关系数 \(R^2=0.9926\)。
如果说这个时候应用多项式插值法，得到的拟合结果是
其相关系数 \(R^2=1\)，我们很好的拟合了该传感器的传感特性，这使得对该传感器的研究更加精确了。
上述只是为了说明插值法的意义，实际上最小二乘法是直线拟合法，不要求所有离散点必须通过该拟合直线；而插值法要求所有离散点必须通过拟合直线。
一般来说，我们将插值问题描述为，
已知 \(n+1\) 个节点 \((x_i,y_i)\)，其中 \(x_i\) 互不相同，不妨设 \(a=x_0<x_1<\cdots<x_n=b\)，那么求任一插值点 \(x^*\neq x_i\) 处的插值 \(y^*\)。
随后构造平面曲线 \(y=\varphi(x)\) 使其通过所有节点，即 \(y_i=\varphi(x_i)\)。
根据上述定义，可以得到插值法的定义为：
设函数 \(y=f(x)\) 在区间 \([a,b]\) 上有定义，且已知在点 \(a=x_0<x_1<\cdots<x_n=b\) 上的值 \(f(x_i)=y_i(i=0,1,\cdots,n)\)，若存在一简单函数 \(\varphi(x)\)，使 \[ \varphi(x_i)=y_i(i=0,1,\cdots,n) \] 成立，就称 \(\varphi(x)\) 为 \(f(x)\) 的插值函数，点 \(x_i(i=0,1,\cdots,n)\) 为插值节点，包括插值节点的区间 \([a,b]\) 称为插值区间，求插值函数 \(\varphi(x)\) 的方法称为插值法。
对于被插函数 \(f(x)\) 和插值函数 \(\varphi(x)\) 在节点 \(x_i\) 处的函数值必然相等，但在节点外 \(\varphi(x)\) 的值可能就会偏离 \(f(x)\)，因此 \(\varphi(x)\) 近似代替 \(f(x)\) 必然存在截断误差。
如果说 \(\varphi_n(x)\) 为次数不超过 \(n\) 的代数多项式 \[ \varphi_n(x)=a_0+a_1x+\cdots+a_nx^n \] 那么就称 \(\varphi_n(x)\) 为插值多项式，相应的插值法称为多项式插值法；
若 \(\varphi_n(x)\) 为分段多项式，就称为分段插值法。
一般来说会选择多项式函数作为插值法的目标插值函数，这是因为其计算简便、性质简单。
多项式插值法
唯一解
由定义可以构造方程组 \[ \begin{cases} a_0+a_1x_0+a_2x_0^2+\cdots+a_nx_0^n=y_0\\ a_0+a_1x_1+a_2x_1^2+\cdots+a_nx_1^n=y_1\\ \vdots\\ a_0+a_1x_n+a_2x_n^2+\cdots+a_nx_n^n=y_n\\ \end{cases} \] 那么该方程组的行列式为 \[ \det{(\mathcal{A})}= \begin{vmatrix} 1&x_0&\cdots&x_0^n\\ 1&x_1&\cdots&x_1^n\\ \vdots&\vdots&\ddots&\vdots\\ 1&x_n&\cdots&x_n^n\\ \end{vmatrix}=\prod(x_i-x_j) \] 只要满足任意 \(x_i\neq x_j\) 且插值多项式 \(\varphi(x)\) 的阶数 \(\leq n\)，那么该方程组就有唯一解。
误差估计
上文提到插值法会造成截断误差，而插值多项式的截断误差可以写为 \[ R_n(x)=f(x)-\varphi_n(x) \] 也称为插值多项式的余项。
截断误差的证明过程是复杂的，这里仅给出结论 \[ R_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}w_{n+1}(x),w_{n+1}(x)=\prod^n_{i=0}(x-x_i) \] 但是这个表达式只有在 \(f(x)\) 的高阶导数存在时才能应用，且 \(\xi\) 无法确定，但可以估算 \(\xi\in (a,b)\)，如果可以求出 \[ \max_{a<x<b}|f^{(n+1)}(x)|=M_{n+1} \] 那么可以得到误差限 \[ |R_n(x)|\leq \frac{M_{n+1}}{(n+1)!}|w_{n+1}(x)| \]
拉格朗日插值法
线性插值
先讨论 \(n=1\) 的情况，\(\varphi_1(x)=a_0+a_1x\)，这也被称为线性插值
那么根据两点，可以构造出方程组 \[ \begin{cases} \varphi_1(x_0)=a_0+a_1x_0=y_0=f(x_0)\\ \varphi_1(x_1)=a_0+a_1x_1=y_1=f(x_1)\\ \end{cases} \] 使用两点式公式，可以得到 \[ \varphi_1(x)=\frac{x-x_1}{x_0-x_1}y_0+\frac{x-x_0}{x_1-x_0}y_1 \] 似乎可以抽象出一个函数 \(l_i(x)\) 将插值多项式表示为 \[ \varphi_1(x)=l_0(x)y_0+l_1(x)y_1 \]
抛物插值
再讨论 \(n=2\) 的情况，\(\varphi_2(x)=a_0+a_1x+a_2x^2\)，这也被称为抛物插值
同样的，可以构造出 \[ \begin{align} \varphi_2(x)&=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}y_0\\ &+\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}y_0\\ &+\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}y_0 \end{align} \] 似乎同样可以抽象出一个函数 \(l_i(x)\) 将插值多项式表示为 \[ \varphi_2(x)=l_0(x)y_0+l_1(x)y_1+l_2(x)y_2 \]
一般形式
更一般地，可以发现插值多项式就是 \(n+1\) 个线性函数的线性组合。
将规律延伸到 \(n\) 次，那么可以得到 \[ L_n(x)=l_0(x)y_0+\cdots+l_n(x)y_n=\sum^n_{i=0}y_il_i(x) \] 其中，\(L_n\) 称为拉格朗日插值函数，\(l_i(x)\) 称为拉格朗日基函数，定义为 \[ l_i(x_j)=\begin{cases} 0,&j\neq i\\ 1,&j=i \end{cases} \] 为了计算方便，可以写出数值形式为 \[ l_i(x)=\frac{(x-x_0)\cdots(x-x_{n})}{(x_i-x_0)\cdots(x_i-x_n)}=\prod^n_{j=0,j\neq i}\frac{x-x_j}{x_i-x_j} \] 而插值多项式为 \[ L_n(x)=\sum^n_{i=0}y_il_i(x)=\sum^n_{i=0}\left(\prod^n_{j=0,j\neq i}\frac{x-x_j}{x_i-x_j}\right) \]
除此之外还有 \(w_{n+1}(x)\) 形式下的拉格朗日插值函数 \[ \begin{align} w_{n+1}(x)&=\prod^n_{j=0}(x-x_j)\\ w_{n+1}'(x_i)&=\prod^n_{j=0,j\neq i}(x_i-x_j)\\ L_n(x)&=\sum^n_{i=0}\frac{w_{n+1}(x)}{(x-x_i)w_{n+1}'(x_i)}y_i \end{align} \]
优缺点
拉格朗日插值法的公式简洁，结构紧凑，理论分析方便，但改变一个数据点时，全部的插值基函数都将发生改变，基函数没有承袭性，且数据点很多时，计算量很大，为 \(2n^2+2n\)。
牛顿插值法
规律延伸
同理首先考虑线性插值 \[ \varphi_1(x)=y_0+\frac{y_1-y_0}{x_1-x_0}(x-x_0) \] 但此时我们使用点斜式，再来观察抛物插值 \[ \varphi_2(x)=y_0+\frac{y_1-y_0}{x_1-x_0}(x-x_0)+\frac{\frac{y_2-y_0}{x_2-x_0}-\frac{y_1-y_0}{x_1-x_0}}{x_2-x_1}(x-x_0)(x-x_1) \] 我们可以从中抽象出一个函数 \(p_i\)，使得 \[ \varphi_n(x)=p_0+p_1(x-x_0)+p_2(x-x_0)(x-x_1)+\cdots+p_n\prod^{n-1}_{i=0}(x-x_i) \] 如果说接着推导 \(p_n\)，那么其形式会十分复杂，故需要引入差商的概念来简化运算。
差商
对函数 \(f(x)\) 有 \(x_0,x_1,\cdots\) 一系列互不相同的点，称 \[ f[x_i,x_j]=\frac{f(x_i)-f(x_j)}{x_i-x_j}(i\neq j) \] 为 \(f(x)\) 关于点 \(x_i, x_j\) 的一阶差商（均差）。
类比斜率。
中括号是为了有意区分函数值和差商的概念，可以认为 \(f[x_i]=f(x_i)\) 是零阶差商。
而称 \[ f[x_i,x_j,x_k]=\frac{f[x_i,x_j]-f[x_j,x_k]}{x_i-x_k}(i\neq j\neq k) \] 为 \(f(x)\) 关于点 \(x_i,x_j,x_k\) 的二阶差商（一阶差商的差商）。
那么 \(k\) 阶差商如下 \[ f[x_0,x_1,\cdots,x_k]=\frac{f[x_0,x_1,\cdots,x_{k-1}]-f[x_1,x_2,\cdots,x_{k}]}{x_0-x_k} \]
差商的性质
差商具有线性性，若 \(f(x)=a\varphi(x)+b\phi(x)\)，则对于任意正整数 \(k\) 都有 \[ f[x_0,x_1,\cdots,x_k]=a\varphi[x_0,x_1,\cdots,x_k]+b\phi[x_0,x_1,\cdots,x_k] \] \(k\) 阶差商可以由函数值的线性组合表示，即 \[ f[x_0,x_1,\cdots,x_k]=\sum^k_{i=0}\frac{f(x_i)}{(x_i-x_0)\cdots(x_i-x_{k})}=\sum^n_{i=0}\frac{f(x_i)}{w_{k+1}'(x_i)} \] 差商具有对称性，即任意调换节点的次序，差商的值不变，如 \[ f[x_0,x_1,x_2]=f[x_0,x_2,x_1]=f[x_2,x_1,x_0] \] 若 \(f(x)\) 是 \(x\) 的 \(n\) 次多项式，则其一阶差商 \(f[x,x_i]\) 是 \(n-1\) 次多项式。
若 \(f(x)\) 在 \([a,b]\) 上存在 \(n\) 阶导数，且 \(x_i\in[a,b]\)，那么 \(n\) 阶差商与导数的关系为 \[ f[x_0,x_1,\cdots,x_n]=\frac{f^{(n)}(\xi)}{n!},\xi\in(a,b) \]
牛顿多项式
用差商的形式表示为 \[ \begin{align} N_n(x)&=f(x_0)\\ &+(x-x0)f[x_0,x_1]\\ &+(x-x_0)(x-x_1)f[x_0,x_1,x_2]\\ &+\cdots\\ &+(x-x_0)\cdots(x-x_{n-1})f[x_0,\cdots,x_n] \end{align} \] 其中牛顿多项式的误差也可以用差商表示 \[ R_n(x)=(x-x_0)\cdots(x-x_n)f[x,x_0,\cdots,x_n]=w_{n+1}(x)f[x,x_0,\cdots,x_n] \] 牛顿插值的优点是，每增加一个节点，插值多项式只增加一项，因此便于递推运算和程序设计，即 \[ N_n(x)=N_{n-1}(x)+f[x_0,\cdots,x_n](x-x_0)\cdots(x-x_{n-1}) \] 且与拉格朗日插值法相比，多项式等价，误差一致。
一般来说，会根据 \(x_i,y_i\) 的值构造差分表，比如说以下数据
| \(x_i\) | \(y_i\) |
|---|---|
| 0.40 | 0.41075 |
| 0.55 | 0.57815 |
| 0.65 | 0.69675 |
| 0.80 | 0.88811 |
| 0.90 | 1.02652 |
| 1.05 | 1.25382 |
我们可以构造这样的差分表
| x_i | y_i | 一阶差商 | 二阶差商 | 三阶差商 | 四阶差商 | 五阶差商 |
|---|---|---|---|---|---|---|
| 0.40 | 0.41075 | |||||
| 0.55 | 0.57815 | 1.11600 | ||||
| 0.65 | 0.69675 | 1.18600 | 0.28000 | |||
| 0.80 | 0.88811 | 1.27573 | 0.35893 | 0.197333 | ||
| 0.90 | 1.02652 | 1.38410 | 0.43348 | 0.212952 | 0.0312381 | |
| 1.05 | 1.25382 | 1.51533 | 0.52493 | 0.228667 | 0.0314286 | 0.000380952 |
其中，对角线即为所求的差商 \(f[x_0],f[x_0,x_1],\cdots,f[x_0,\cdots,x_n]\)。
差分
定义
向前差分，Forward Difference（F-差分） \[ \Delta f_k=f_{k+1}-f_k \] 向后差分，Backward Difference（B-差分） \[ \nabla f_k=f_k-f_{k-1} \]
性质
对于向前差分而言，具有连续性 \[
\begin{align}
\Delta f_k&=f_{k+1}-f_k\\
\Delta^2 f_k&=\Delta f_{k+1}-\Delta f_k\\
&\vdots\\
\Delta^m f_k&=\Delta^{m-1} f_{k+1}-\Delta^{m-1} f_k
\end{align}
\] 
线性性 \[ f(x)=a\varphi(x)+b\phi(x)\to \Delta^m f_k=a\Delta^m \varphi_k+b\Delta^m \phi_k \] 线性表示 \[ \Delta^m f_k=\sum^m_{j=0}(-1)^j\binom{j}{m} f_{k+m-j} \] 差分与差商的关系 \[ f[x_k,\cdots,x_{k+m}]=\frac{1}{m!}\frac{1}{h^m}\Delta^m f_k \] 差分与导数的关系 \[ f^{(m)}(\xi)=\frac{\Delta^m f_k}{h^m} \]
等距节点牛顿插值法
实际使用中，通常我们得到的是等距节点下的数据值，例如说 \((0,1.002),(1,1.208),\cdots,(x_n,y_n)\)，满足 \(h=x_i-x_{i-1}\) 为常数（等差）。
不妨将差分与差商的关系代入牛顿插值公式，可以得到 \[ N_n(x)=f_0+\frac{\delta f_0}{h}w_1(x)+\frac{\Delta^2 f_0}{2!h^2}w_2(x)+\cdots+\frac{\Delta^n f_0}{n!h^n} \] 其中，\(h\) 为等距节点的等差。
为方便起见，将 \(w_{n+1}(x)\) 化为无关 \(x\) 的函数，即不妨取 \(x=x_0+th\)，其中 \(t\) 为间距倍数，使得 \[ w_{n+1}(x)=w_{n+1}(x_0+th)=(t-0)(t-1)\cdots(t-n)h^{n+1} \] 那么牛顿多项式化为 F-插值公式 \[ N_n(x)=f_0+t\Delta f_0+\frac{t(t-1)}{2!}\Delta^2 f_0+\cdots+\frac{t(t-1)\cdots(t-n+1)}{n!}\Delta^n f_0 \] 其余项（误差）化为 \[ R_n(x)=\frac{t(t-1)\cdots(t-n)}{(n+1)!}h^{n+1}f^{(n+1)}(\xi) \] 误差限为 \[ |R_n(x)|\leq \left|\frac{t(t-1)\cdots(t-n)}{(n+1)!}h^{n+1}M_{n+1}\right|,|f^{(n+1)}(\xi)|\leq M_{n+1} \] 同理，B-差分也有类似的关系，可以得到 B-插值公式 \[ N_n(x)=f_n+t\nabla f_n+\frac{t(t-1)}{2!}\nabla^2 f_n+\cdots+\frac{t(t-1)\cdots(t-n+1)}{n!}\nabla^n f_n \] 当插值点 \(x\) 接近数据表头时，一般采用向前插值法；
当插值点 \(x\) 接近数据表尾时，一般采用向后插值法。
看的是 \(x\) 接近 \(f_0\) 还是 \(f_n\)，越相近截断误差越小。
埃尔米特插值法
埃尔米特插值多项式
许多实际的插值问题不但要求节点处函数值相同，还要求它的函数有相同的一阶、二阶导数。
满足这种要求的插值多项式就是埃尔米特插值多项式，其插值多项式有 \(2n+1\) 次数，即 \[ H_{2n+1}(x)=a_0+a_1x+\cdots+a_{2n+1}x^{2n+1} \] 故需要 \(2n+2\) 个条件（数据点），其满足 \[ H_{2n+1}(x_i)=f(x_i),H'_{2n+1}(x_i)=f'(x_i) \] 余项为 \[ R(x)=\frac{f^{2n+2}(\xi)}{(2n+2)!}w_{n+1}^2(x) \]