多项式插值

多项式插值

1、多项式插值的定义

给定$n+1$个点$\{(x_i,y_i)\}^{n}_{i=0}$（称为 插值点），所谓 多项式插值 就是找到一个多项式（称为 插值多项式）

$$y=P(x)=a_kx^k+a_{k-1}x^{k-1}+\cdots+a_1x+a_0$$
使得它满足条件
$$y_i=P(x_i)，其中i=0,1,\dots,n$$
也就是说，多项式 $y=P(x)$的图像要经过给定的 $n+1$个点。

在实际应用中，这些插值点可能来自某次实验测量所得的数据，也可能来自某个复杂函数$y=f(x)$的值。通过计算插值多项式，我们可以找到这些实验数据间的规律，或者使用简单的多项式函数$y=P(x)$来近似复杂的函数$y=f(x)$。

2、插值多项式的存在唯一性和误差

满足以上条件的多项式总是存在的。特别地，我们可以证明：

定理一：给定$n+1$个点$\{(x_i,y_i)\}^{n}_{i=0}$，若$x_i$两两不同，则存在 唯一 一个次数不超过$n$的多项式$y=P(x)$，使得$y_i=P(x_i)$（$i=0,1,\dots,n$）成立。
证明：利用范德蒙德矩阵和代数学基本定理即得。

注意：在本文中我们均假定$x_i$互不相同。

当$y_i$的值来自某个函数$y=f(x)$，且$f(x)$具有$n+1$阶连续导数时，下面的定理可以用来计算多项式插值的（截断）误差。

定理二：给定$n+1$个点$\{(x_i,y_i)\}^{n}_{i=0}$，其中$y_i=f(x_i)$，进一步假设函数$f(x)$具有$n+1$阶连续导数，则插值多项式$P(x)$的误差$R(x)$为

$$R(x)=f(x)-P(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}(x-x_0)(x-x_1)\cdots(x-x_n)$$
其中 $\min\{x_i\}\leq\xi\leq\max\{x_i\}$。

3、插值多项式的计算

给定$n+1$个点$\{(x_i,y_i)\}^{n}_{i=0}$，计算插值多项式的主要方法有：直接法、拉格朗日多项式插值和牛顿多项式插值。下面我们分别介绍这三种方法。（注意，根据定理一，这三种方法得到的插值多项式在理论上说应该是一致的，而且误差也相同。）

3.1 直接法

根据定理一，假设插值多项式为

$$y=P(x)=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0$$
由插值条件 $y_i=P(x_i)$（ $i=0,1,\dots,n$)，我们得到关于系数 $a_n,a_{n-1},\dots,a_1,a_0$的线性方程组
$$\left(\begin{array}{lllll} x^n_0 & x^{n-1}_0 & \cdots & x_0 & 1\\ x^n_1 & x^{n-1}_1 & \cdots & x_1 & 1\\ \vdots & \vdots & \ddots & \vdots & \vdots \\ x^n_{n-1} & x^{n-1}_{n-1} & \cdots & x_{n-1} & 1\\ x^n_n & x^{n-1}_n & \cdots & x_n & 1\\ \end{array}\right) \left(\begin{array}{l} a_n \\ a_{n-1}\\ \vdots \\ a_1 \\ a_0\\ \end{array}\right)=\left(\begin{array}{l} y_0 \\ y_1\\ \vdots \\ y_{n-1} \\ y_n\\ \end{array}\right)$$
通过求解这个线性方程组，即得到插值多项式。

优点：直接，性质一目了然。
缺点：待求解的线性方程组的系数矩阵为 范德蒙德 (Vandermonde)矩阵，它是一个病态矩阵，这使得在实际求解方程组时将产生很大的误差。

3.2 拉格朗日多项式插值

拉格朗日(Lagrange)多项式插值的原理是：先构造一组 拉格朗日基函数 $L_i(x)$（$i=0,1,\dots,n$），这些基函数为次数不超过$n$的多项式，且具有性质

$$L_i(x_j)=\delta_{ij}=\begin{cases}1 & \mbox{当 }i=j\\0 & \mbox{当 }i\neq j\end{cases}$$
然后将这些基函数做线性组合，得到拉格朗日插值多项式

$$L(x)=\sum^n_{i=0}L_i(x)y_i$$
容易验证，多项式 $L(x)$满足插值条件 $y_i=L(x_i)$（ $i=0,1,\dots,n$)。

拉格朗日基函数$L_i(x)$的构造如下：

由基函数的性质，当$j\in\{0,1,\dots,n\}\setminus\{i\}$时$L_i(x_j)=0$，即$x_j$为$L_i(x)$的零点，可以假设

$$L_i(x)=K\cdot \prod_{j\neq i}(x-x_j)$$
其中， $K$为待定系数。再由$L_i(x_i)=1$，得到
$$L_i(x_i)=K\cdot \prod_{j\neq i}(x_i-x_j)=1$$
从而得到
$$K=\frac{1}{\prod_{j\neq i}(x_i-x_j)}$$
因此，基函数
$$L_i(x)=\frac{\prod_{j\neq i}(x-x_j)}{\prod_{j\neq i}(x_i-x_j)}$$
令 $\omega(x)=(x-x_0)(x-x_1)\cdots(x-x_n)$，则 $L_i(x)$还可以表示为
$$L_i(x)=\frac{\omega_{n+1}(x)}{(x-x_i)\omega'_{n+1}(x_i)}$$

下面的定理说明$L_i(x)$称为基函数的原因。

定理三：令$P_n(x)$为全体次数不超过$n$的多项式构成的集合，则$\{L_i(x)\}^n_{i=0}$是线性空间$P_n(x)$的一组基。

Matlab代码

function [y,Lb] = LagrangeInterpolation(X,Y,x)
% 拉格朗日多项式插值函数
% 注意：插值点的个数为n，差值多项式的次数为n-1
%
% 输入参数
% X,Y: 插值点坐标
% x: 求值点
%
% 输出参数
% y: 拉格朗日插值多项式在x点的值
% Lb: 拉格朗日基函数在x点的值

if length(X) ~= length(Y)
    error('X和Y的长度不相等');
end
n = length(X); %获取插值点的个数
%初始化
y = 0;
Lb = ones(1,n);
for i = 1:n
    for j = 1:n %计算拉格朗日基函数在x点的值
        if j ~= i
            Lb(i) = Lb(i) * (x-X(j))/(X(i)-X(j));
        end
    end
    y = y + Lb(i)*Y(i); %计算拉格朗日插值多项式的值
end

end

3.3 均差与牛顿多项式插值

牛顿多项式插值是基于均差的计算。首先定义均差如下。

函数$f(x)$关于点$x_i,x_{i+1}$的 一阶均差（或差商）为

$$f[x_i,x_{i+1}]=\frac{f(x_{i+1})-f(x_{i})}{x_{i+1}-x_i}$$
一阶均差 $f[x_i,x_{i+1}]$反映了函数 $y=f(x)$在区间 $[x_i,x_{i+1}]$的平均变化率。用递归的方式，我们定义 二阶均差 为
$$f[x_i,x_{i+1},x_{i+2}]=\frac{f[x_{i+1},x_{i+2}]-f[x_{i},x_{i+1}]}{x_{i+2}-x_i}$$
同理， $k$阶均差 为
$$f[x_i,x_{i+1},\dots,x_{i+k}]=\frac{f[x_{i+1},x_{i+2},\dots,x_{i+k}]-f[x_{i},x_{i+1},\dots,x_{i+k-1}]}{x_{i+k}-x_i}$$
特别地，$0$阶均差 定义为$f[x_i]=f(x_i)$。

根据均差的定义，构造均差表如下。

$x_k$	$0$阶	$1$阶	$2$阶	$3$阶	$\cdots$	因子
$x_0$	$f(x_0)$					$1$
$x_1$	$f(x_1)$	$f[x_0,x_1]$				$x-x_0$
$x_2$	$f(x_2)$	$f[x_1,x_2]$	$f[x_0,x_1,x_2]$			$(x-x_0)(x-x_1)$
$x_3$	$f(x_3)$	$f[x_2,x_3]$	$f[x_1,x_2,x_3]$	$f[x_0,x_1,x_2,x_3]$		$\prod^2_{i=0}(x-x_i)$
$\vdots$	$\vdots$	$\vdots$	$\vdots$	$\vdots$	$\ddots$	$\vdots$

如果将$x$也看作一个点，由均差的定义可以得到

$$\begin{eqnarray*} f(x)&=&f(x_0)+f[x,x_0](x-x_0)\\ f[x,x_0]&=&f[x_0,x_1]+f[x,x_0,x_1](x-x_1)\\ f[x,x_0,x_1]&=&f[x_0,x_1,x_2]+f[x,x_0,x_1,x_2](x-x_2)\\ &\vdots&\\ f[x,x_0,x_1,\dots,x_{n-1}]&=&f[x_0,x_1,\dots,x_n]+f[x,x_0,x_1,\dots,x_n](x-x_n) \end{eqnarray*}$$
把以上各式由下而上逐步代入，得到

$$\begin{eqnarray*} f(x)&=&f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)\\ &&+\cdots+f[x_0,x_1,\dots,x_n](x-x_0)(x-x_1)\cdots(x-x_{n-1})\\ &&+f[x,x_0,x_1,\dots,x_n](x-x_0)(x-x_1)\cdots(x-x_n) \end{eqnarray*}$$
其中，
$$\begin{eqnarray*} N(x)&=&f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)\\ &&+\cdots+f[x_0,x_1,\dots,x_n](x-x_0)(x-x_1)\cdots(x-x_{n-1}) \end{eqnarray*}$$
称为牛顿插值多项式。
$$R(x)=f[x,x_0,x_1,\dots,x_n](x-x_0)(x-x_1)\cdots(x-x_n)$$
为插值余项。

由 定理一 和 定理二 得到均差和导数的关系如下

$$f[x,x_0,x_1,\dots,x_n]=\frac{f^{(n+1)}(\xi)}{(n+1)!}$$
其中 $\min\{x_i\}\leq\xi\leq\max\{x_i\}$。

Matlab代码

function [y,Nt]=NewtonInterpolation(X,Y,x)
% 牛顿多项式插值函数
% 注意：插值点的个数为n，差值多项式的次数为n-1
%
% 输入参数
% X,Y: 插值点坐标
% x: 求值点
%
% 输出参数
% y：牛顿差值多项式在x点的值
% Nt：均差表

if length(X) ~= length(Y)
    error('X和Y的长度不相等');
end
n = length(X); 
Nt = zeros(n); %初始化均差表，按列存放各阶均差
Nt(1,1) = Y(1); %0阶均差
for i = 2:n %按行计算均差表
    Nt(i,1) = Y(i); %0阶均差
    for j = 2:i
        Nt(i,j) = (Nt(i,j-1)-Nt(i-1,j-1))/(X(i)-X(i-j+1));
    end
end
%计算牛顿插值多项式在x点上的值
w = 1;
y = Nt(1,1) * w;
for i = 2 : n
    w = w * (x - X(i-1));
    y = y + Nt(i,i) * w;
end

end

3.4 拉格朗日多项式插值和牛顿多项式插值的比较

拉格朗日多项式插值的计算量大于牛顿多项式插值的计算量。特别地，当新增一个插值点时，拉格朗日插值需要重新计算全部的基函数，而牛顿插值只需计算均差表中新的一行的值即可。