数值分析——牛顿插值多项式
数值分析——牛顿插值多项式
牛顿插值多项式是一种逐次生成的插值函数，特别适用于需要频繁增减插值点的场景。本文将详细介绍牛顿插值多项式的理论推导、均差的概念以及牛顿前插公式的推导，并提供MATLAB实现代码及测试用例。
理论推导
1.插值多项式的逐次生成
基于拉格朗日插值的推导思路，由一阶的线性插值函数来逐次深入到高阶的牛顿插值多项式。首先需要了解到线性插值的另一种表现形式——点斜式：
$$
P_1(x) = L_1(x) = y_0 + \frac{y_1 - y_0}{x_1 - x_0}(x - x_0) = P_0(x) + a_1(x - x_0) \tag{1}
$$
其中$P_0(x) = y_0$，且$a_1 = \frac{y_1 - y_0}{x_1 - x_0} = \frac{f(x_1) - f(x_0)}{x_1 - x_0}$被定义为函数的差商。点斜式的插值函数也满足：
$$
P_1(x_0) = y_0, \quad P_1(x_1) = y_1 \tag{2}
$$
考虑三个点的抛物线插值函数$L_2(x)$满足：
$$
P_2(x_0) = y_0, \quad P_2(x_1) = y_1, \quad P_2(x_2) = y_2 \tag{3}
$$
基于此条件可以发现线性插值函数$P_1(x)$，其实就已经满足公式(3)的前两个条件，只需要在其基础上增加一个多项式使$P_2(x)$满足公式(3)的第三个条件即可：
$$
P_2(x) = P_1(x) + a_2(x - x_0)(x - x_1) \tag{4}
$$
其中$a_2$为待定系数可通过$P_2(x_2) = y_2$求得：
$$
a_2 = \frac{P_2(x_2) - P_1(x_2)}{(x - x_0)(x - x_1)} = \frac{\frac{f(x_2) - f(x_1)}{x_2 - x_1} - \frac{f(x_1) - f(x_0)}{x_1 - x_0}}{x_2 - x_1} \tag{5}
$$
2.均差
可以看到系数的表达式十分复杂，为了推导与理解，这里引入均差（差商）的定义：
$$
f[x_0, x_1, \cdots, x_{k-1}, x_k] = \frac{f[x_0, x_1, \cdots, x_{k-2}, x_k] - f[x_0, x_1, \cdots, x_{k-2}, x_{k-1}]}{x_k - x_{k-1}} = \frac{f[x_1, \cdots, x_k] - f[x_0, \cdots, x_{k-1}]}{x_k - x_0} \tag{6}
$$
因此$f[x_0, x_k] = \frac{f(x_k) - f(x_0)}{x_k - x_0}$被称为函数$f(x)$关于点$x_0$、$x_k$的一阶均差，$f[x_0, x_1, x_k] = \frac{f[x_0, x_k] - f[x_0, x_1]}{x_k - x_1}$被称为函数$f(x)$关于点$x_0$、$x_1$、$x_k$的二阶均差，公式(6)被称为函数$f(x)$关于点$x_0$、$x_1$、$\cdots$、$x_{k-1}$、$x_k$的$k$阶均差。公式(6)的另外一个表达形式为：
$$
f[x_0, x_1, \cdots, x_{k-1}, x_k] = \sum_{j=0}^k \frac{f(x_j)}{(x_j - x_0) \cdots (x_j - x_{j-1})(x_j - x_{j+1}) \cdots (x_j - x_k)} \tag{7}
$$
根据文献，公式(6)到公式(7)的推导可以通过归纳法证明。对于在区间内存在$n$阶导数的函数，$n$阶均差可以表示为：
$$
f[x_0, x_1, \cdots, x_{k-1}, x_k] = \frac{f^{(n)}(\xi)}{n!}, \quad \xi \in [a, b] \tag{8}
$$
3.牛顿插值多项式
引入均差的定义后线性插值函数表达式就可以推导为：
$$
P_1(x) = f(x_0) + f[x_0, x_1](x - x_0) \tag{9}
$$
抛物线插值函数表达式就可以推导为：
$$
P_2(x) = P_1(x) + f[x_0, x_1, x_2](x - x_0)(x - x_1) \tag{10}
$$
$n$阶插值函数可以表示为：
$$
\begin{align}
P_n(x) &= P_{n-1}(x) + f[x_0, x_1, \cdots, x_{n-1}, x_n](x - x_0)(x - x_1) \cdots (x - x_{n-1}) \
&= P_{n-2}(x) + f[x_0, x_1, \cdots, x_{n-2}, x_{n-1}](x - x_0)(x - x_1) \cdots (x - x_{n-2}) \
&\quad + f[x_0, x_1, \cdots, x_{n-1}, x_n](x - x_0)(x - x_1) \cdots (x - x_{n-1}) \
&= \cdots \
&= f(x_0) + f[x_0, x_1](x - x_0) + \cdots + f[x_0, x_1, \cdots, x_{n-1}, x_n](x - x_0)(x - x_1) \cdots (x - x_{n-1}) \tag{11}
\end{align}
$$
该式就被称为牛顿插值多项式。
4.牛顿前插公式
在上文中的情况中，插值点都是随机分布。在一些特定的场合下，插值点的间隔是固定的，此时牛顿插值多项式就可以转化为即将介绍的牛顿前插公式。
由于插值点之间是等距的，故设：
$$
x_k = x_0 + kh \quad k=0,1,\cdots,n-1,n \tag{12}
$$
其中$h$、$x_0$分别是步长和初始点。设$x_k$点的函数值为$f_k = f(x_k) (k=0,1,\cdots,n-1,n)$，称$\Delta f_k = f_{k+1} - f_{k}$为$x_k$处以$h$为步长的一阶前向差分，称$\Delta^2 f_k = \Delta f_{k+1} - \Delta f_{k}$为$x_k$处以$h$为步长的二阶前向差分，进而
$$
\Delta^n f_k = \Delta^{n-1} f_{k+1} - \Delta^{n-1} f_{k} \tag{13}
$$
称为$x_k$处以$h$为步长的$n$阶前向差分，各阶差分的计算路径如图：
导出均差与差分的关系：
$$
f[x_k, \cdots, x_{k+m}] = \frac{1}{m!} \frac{1}{h^m} \Delta^m f_k, \quad m=1,2,\cdots,n \tag{14}
$$
由此可得牛顿前插公式：
$$
P_n(x_0 + th) = 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
$$
MATLAB实现
1.利用暴力递归法实现牛顿插值多项式
由于牛顿多项式的计算每一步都可以使用前一步的计算结果，可以使用循环或递归的方式进行编程：
% 牛顿多项式插值
function [aArr] = func_NewtonInterpolation(xArr, yArr, n)
% 识别误差
xLength = length(xArr);
yLength = length(yArr);
% 插值数组不匹配
if xLength ~= yLength
error('Error: The length of xArr=%d and yArr=%d are not equal in length.', xLength, yLength);
elseif xLength ~= n
error('Error: The length of xArr=%d and n=%d are not equal.', xLength, n);
elseif yLength ~= n
error('Error: The length of yArr=%d and n=%d are not equal.', yLength, n);
elseif n == 0
error('Error: The n should gainer than %d.', n);
end
% 递归计算牛顿插值多项式
% 底层返回
if n == 1
aArr = yArr(n);
return;
end
% 计算新增多项式
syms x;
newPoly = func_DifferenceQuotient1(xArr, yArr, n);
for i = 1:1:n-1
newPoly = newPoly * (x - xArr(i));
end
% 递归
aArr = sym2poly(poly2sym(func_NewtonInterpolation(xArr(1:1:n-1), yArr(1:1:n-1), n-1), x) + newPoly);
end
2.牛顿插值多项式测试
针对已有函数表：
设计测试程序：
% 测试牛顿插值多项式函数——func_NewtonInterpolation
clc;
clear;
close all;
xArr = [0.40, 0.55, 0.65, 0.80, 0.90, 1.05];
yArr = [0.41075, 0.57815, 0.69675, 0.88811, 1.02652, 1.25382];
polyval(func_NewtonInterpolation(xArr(1:1:5), yArr(1:1:5), 5), 0.596)
测试结果：
3.利用暴力递归法实现牛顿前插多项式
实现代码如下：
% 差分形式的牛顿多项式插值
function [aArr] = func_DifferenceNewtonInterpolation(yArr, n)
% 识别误差
yLength = length(yArr);
if yLength ~= n
error('Error: The length of yArr=%d and n=%d are not equal.', yLength, n);
elseif n == 0
error('Error: The n should gainer than %d.', n);
end
% 递归计算牛顿插值多项式
% 底层返回
if n == 1
aArr = yArr;
return;
end
% 计算新增多项式
syms t;
newpoly = func_Difference(yArr, n, n-1) / factorial(n-1);
for i = 1:1:n-1
newpoly = newpoly * (t - i + 1);
end
% 递归
aArr = sym2poly(poly2sym(func_DifferenceNewtonInterpolation(yArr(1:1:n-1), n-1), t) + newpoly);
end
4.牛顿前插多项式测试
针对$\cos(x)$已有函数表：
设计测试程序：
% 测试前插牛顿插值函数——func_DifferenceNewtonInterpolation
clc;
clear;
close all;
xArr = 0.00:0.10:0.50;
yArr = [1, 0.995, 0.98007, 0.95534, 0.92106, 0.87758];
polyval(func_DifferenceNewtonInterpolation(yArr(1:1:5), 5), 0.048)
cos(0.048)
测试结果：
总结
以上在本文中对牛顿插值多项式与牛顿前插多项式的推导进行了简单的介绍，并提供了MATLAB的实现代码与测试用例。
参考文献
[1] 李庆扬, 王能超, 易大义. 数值分析. 第5版 [M]. 清华大学出版社, 2008.