数值分析——牛顿插值多项式

关键字：Matalb；逐次生成；曲线拟合

系列文章目录

数值分析——拉格朗日插值
数值分析——埃尔米特（Hermit）插值
数值分析——分段低次插值
数值分析——三次样条插值

---

---

前言

  在之前的文章中简单介绍了线性插值、抛物线插值、拉格朗日插值，以上的插值方法推导简单、公式结构紧凑是许多曲线拟合的首选方式。但是如果遇到增减插值点的情况时，需要全部进行重新计算。为了计算的方便，逐次生成插值多项式的插值函数——牛顿插值多项式便应运而生。

本文对一些复杂的流程与概念进行了简化与省略，如果想详细了解的小伙伴可以参考李庆扬老师的《数值分析 -第5版》的第2章（需要电子版请私戳）

---

一、理论推导

1.插值多项式的逐次生成

  基于拉格朗日插值的推导思路，由一阶的线性插值函数来逐次深入到高阶的牛顿插值多项式。首先需要了解到线性插值的另一种表现形式——点斜式：
$$P_1\left(x\right)=L_1\left(x\right)=y_0+\frac{y_1-y_0}{x_1-x_0}\left(x-x_0\right)=P_0\left(x\right)+a_1\left(x-x_0\right)\text{\hspace{1em}(1)}$$其中$P_0\left(x\right)=y_0$，且$a_1=\frac{y_1-y_0}{x_1-x_0}=\frac{f\left(x_1\right)-f\left(x_0\right)}{x_1-x_0}$被定义为函数的差商。点斜式的插值函数也满足：
$$P_1\left(x_0\right)=y_0,P_1\left(x_1\right)=y_1\text{\hspace{1em}(2)}$$考虑三个点的抛物线插值函数$L_2\left(x\right)$满足：
$$P_2\left(x_0\right)=y_0,P_2\left(x_1\right)=y_1,P_2\left(x_2\right)=y_2\text{\hspace{1em}(3)}$$基于此条件可以发现线性插值函数$P_1\left(x\right)$，其实就已经满足公式(3)的前两个条件，只需要在其基础上增加一个多项式使$P_2\left(x\right)$满足公式(3)的第三个条件即可：
$$P_2\left(x\right)=P_1\left(x\right)+a_2\left(x-x_0\right)\left(x-x_1\right)\text{\hspace{1em}(4)}$$其中$a_2$为待定系数可通过$P_2\left(x_2\right)=y_2$求得 ：
$$a_2=\frac{P_2\left(x_2\right)-P_1\left(x_2\right)}{\left(x-x_0\right)\left(x-x_1\right)}=\frac{\frac{f\left(x_2\right)-f\left(x_1\right)}{x_2-x_1}-\frac{f\left(x_1\right)-f\left(x_0\right)}{x_1-x_0}}{x_2-x_1}\text{\hspace{1em}(5)}$$

2.均差

  可以看到系数的表达式十分复杂，为了推导与理解，这里引入均差（差商）的定义：
$$\begin{array}{rl}f\left[x_0,x_1,\cdots ,x_{k-1},x_k\right]& =\frac{f\left[x_0,x_1,\cdots ,x_{k-2},x_k\right]-f\left[x_0,x_1,\cdots ,x_{k-2},x_{k-1}\right]}{x_k-x_{k-1}}\\ & =\frac{f\left[x_1,\cdots ,x_k\right]-f\left[x_0,\cdots ,x_{k-1}\right]}{x_k-x_0}\end{array}\text{\hspace{1em}(6)}$$因此$f\left[x_0,x_k\right]=\frac{f\left(x_k\right)-f\left(x_0\right)}{x_k-x_0}$被称为函数$f\left(x\right)$关于点$x_0$、$x_k$的一阶均差，$f\left[x_0,x_1,x_k\right]=\frac{f\left[x_0,x_k\right]-f\left[x_0,x_1\right]}{x_k-x_1}$被称为函数$f\left(x\right)$关于点$x_0$、$x_1$、$x_k$的二阶均差，公式(6)被称为函数$f\left(x\right)$关于点$x_0、x_1、\cdots 、x_{k-1}、x_k$的k阶均差。公式(6)的另外一个表达形式为：
$$f\left[x_0,x_1,\cdots ,x_{k-1},x_k\right]=\sum _{j=0}^k\frac{f\left(x_j\right)}{\left(x_j-x_0\right)\cdots \left(x_j-x_{j-1}\right)\left(x_j-x_{j+1}\right)\cdots \left(x_j-x_k\right)}\text{\hspace{1em}(7)}$$根据文献，公式(6)到公式(7)的推导可以通过归纳法证明，当然不是不会😉 这里博主就不赘述了。对于在区间内存在$n$阶导数的函数，$n$阶均差可以表示为：
$$f\left[x_0,x_1,\cdots ,x_{k-1},x_k\right]=\frac{f^{\left(n\right)}\left(\xi \right)}{n!},\xi \in \left[a,b\right]\text{\hspace{1em}(8)}$$

3.牛顿插值多项式

  引入均差的定义后线性插值函数表达式就可以推导为：
$$P_1\left(x\right)=f\left(x_0\right)+f\left[x_0,x_1\right]\left(x-x_0\right)\text{\hspace{1em}(9)}$$抛物线插值函数表达式就可以推导为：
$$P_2\left(x\right)=P_1\left(x\right)+f\left[x_0,x_1,x_2\right]\left(x-x_0\right)\left(x-x_1\right)\text{\hspace{1em}(10)}$$$n$阶插值函数可以表示为：
$$\begin{array}{rl}{P}_n\left(x\right)& =P_{n-1}\left(x\right)+f\left[x_0,x_1,\cdots ,x_{n-1},x_n\right]\left(x-x_0\right)\left(x-x_1\right)\cdots \left(x-x_{n-1}\right)\\ & =P_{n-2}\left(x\right)+f\left[x_0,x_1,\cdots ,x_{n-2},x_{n-1}\right]\left(x-x_0\right)\left(x-x_1\right)\cdots \left(x-x_{n-2}\right)+f\left[x_0,x_1,\cdots ,x_{n-1},x_n\right]\left(x-x_0\right)\left(x-x_1\right)\cdots \left(x-x_{n-1}\right)\\ & =\cdots \\ & =f\left(x_0\right)+f\left[x_0,x_1\right]\left(x-x_0\right)+\cdots +f\left[x_0,x_1,\cdots ,x_{n-1},x_n\right]\left(x-x_0\right)\left(x-x_1\right)\cdots \left(x-x_{n-1}\right)\end{array}\text{\hspace{1em}(11)}$$该式就被称为牛顿插值多项式。

4.牛顿前插公式

  在上文中的情况中，插值点都是随机分布。在一些特定的场合下，插值点的间隔是固定的，此时牛顿插值多项式就可以转化为即将介绍的牛顿前插公式。
  由于插值点之间是等距的，故设：
$$x_k=x_0+khk=0,1,\cdots ,n-1,n\text{\hspace{1em}(12)}$$其中$h$、$x_0$分别是步长和初始点。设$x_k$点的函数值为$f_k=f\left(x_k\right)\left(k=0,1,\cdots ,n-1,n\right)$，称$\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\text{\hspace{1em}(13)}$$称为$x_k$处以$h$为步长的$n$阶前向差分，各阶差分的计算路径如图：

图1.1 差分计算路径

导出均差与差分的关系：
$$f\left[x_k,\cdots ,x_{k+m}\right]=\frac{1}{m!}\frac{1}{h^m}{\Delta }^m{f}_k,m=1,2,\cdots ,n\text{\hspace{1em}(14)}$$由此可得牛顿前插公式：
$$P_n\left(x_0+th\right)=f_0+t\Delta f_0+\frac{t\left(t-1\right)}{2!}{\Delta }^2{f}_0+\cdots +\frac{t\left(t-1\right)\cdots \left(t-n+1\right)}{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.牛顿插值多项式测试

  针对已有函数表：

图2.1 函数表

设计测试程序：

% 测试牛顿插值多项式函数——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)

测试结果：

图2.2 牛顿插值多项式计算结果

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\left(x\right)$已有函数表：

图2.3 函数表

设计测试程序：

% 测试前插牛顿插值函数——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)

测试结果：

图2.4 牛顿前插多项式计算结果

---

总结

  以上在本文中对牛顿插值多项式与牛顿前插多项式的推导进行了简单的介绍，并提供了MATLAB的实现代码与测试用例。

参考文献

[1]李庆扬,王能超,易大义. 数值分析.第5版[M]. 清华大学出版社, 2008.