多项式插值

1、多项式插值的定义

给定n+1个点{(xi,yi)}ni=0（称为 插值点），所谓 多项式插值 就是找到一个多项式（称为 插值多项式）

y=P(x)=akxk+ak−1xk−1+⋯+a1x+a0

使得它满足条件

yi=P(xi)，其中i=0,1,…,n

也就是说，多项式y=P(x)的图像要经过给定的n+1个点。

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

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

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

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

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

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

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

R(x)=f(x)−P(x)=f(n+1)(ξ)(n+1)!(x−x0)(x−x1)⋯(x−xn)

其中min{xi}≤ξ≤max{xi}。

3、插值多项式的计算

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

3.1 直接法

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

y=P(x)=anxn+an−1xn−1+⋯+a1x+a0

由插值条件yi=P(xi)（i=0,1,…,n)，我们得到关于系数an,an−1,…,a1,a0的线性方程组

⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜xn0xn1⋮xnn−1xnnxn−10xn−11⋮xn−1n−1xn−1n⋯⋯⋱⋯⋯x0x1⋮xn−1xn11⋮11⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎛⎝⎜⎜⎜⎜⎜⎜⎜anan−1⋮a1a0⎞⎠⎟⎟⎟⎟⎟⎟⎟=⎛⎝⎜⎜⎜⎜⎜⎜⎜y0y1⋮yn−1yn⎞⎠⎟⎟⎟⎟⎟⎟⎟

通过求解这个线性方程组，即得到插值多项式。

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

3.2 拉格朗日多项式插值

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

Li(xj)=δij={10当 i=j当 i≠j

然后将这些基函数做线性组合，得到拉格朗日插值多项式

L(x)=∑i=0nLi(x)yi

容易验证，多项式L(x)满足插值条件yi=L(xi)（i=0,1,…,n)。

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

由基函数的性质，当j∈{0,1,…,n}∖{i}时Li(xj)=0，即xj为Li(x)的零点，可以假设

Li(x)=K⋅∏j≠i(x−xj)

其中，K为待定系数。再由Li(xi)=1，得到

Li(xi)=K⋅∏j≠i(xi−xj)=1

从而得到

K=1∏j≠i(xi−xj)

因此，基函数

Li(x)=∏j≠i(x−xj)∏j≠i(xi−xj)

令ω(x)=(x−x0)(x−x1)⋯(x−xn)，则Li(x)还可以表示为

Li(x)=ωn+1(x)(x−xi)ω′n+1(xi)

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

定理三：令Pn(x)为全体次数不超过n的多项式构成的集合，则{Li(x)}ni=0是线性空间Pn(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的长度不相等');endn = 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); %计算拉格朗日插值多项式的值endend

3.3 均差与牛顿多项式插值

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

函数f(x)关于点xi,xi+1的 一阶均差（或差商）为

f[xi,xi+1]=f(xi+1)−f(xi)xi+1−xi

一阶均差f[xi,xi+1]反映了函数y=f(x)在区间[xi,xi+1]的平均变化率。用递归的方式，我们定义 二阶均差 为

f[xi,xi+1,xi+2]=f[xi+1,xi+2]−f[xi,xi+1]xi+2−xi

同理，k阶均差 为

f[xi,xi+1,…,xi+k]=f[xi+1,xi+2,…,xi+k]−f[xi,xi+1,…,xi+k−1]xi+k−xi

特别地，0阶均差 定义为f[xi]=f(xi)。

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

xk 0阶 1阶 2阶 3阶 ⋯ 因子 x0 f(x0) 1 x1 f(x1) f[x0,x1] x−x0 x2 f(x2) f[x1,x2] f[x0,x1,x2] (x−x0)(x−x1) x3 f(x3) f[x2,x3] f[x1,x2,x3] f[x0,x1,x2,x3] ∏2i=0(x−xi) ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ ⋮

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

f(x)f[x,x0]f[x,x0,x1]f[x,x0,x1,…,xn−1]===⋮=f(x0)+f[x,x0](x−x0)f[x0,x1]+f[x,x0,x1](x−x1)f[x0,x1,x2]+f[x,x0,x1,x2](x−x2)f[x0,x1,…,xn]+f[x,x0,x1,…,xn](x−xn)

把以上各式由下而上逐步代入，得到

f(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+⋯+f[x0,x1,…,xn](x−x0)(x−x1)⋯(x−xn−1)+f[x,x0,x1,…,xn](x−x0)(x−x1)⋯(x−xn)

其中，

N(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+⋯+f[x0,x1,…,xn](x−x0)(x−x1)⋯(x−xn−1)

称为牛顿插值多项式。

R(x)=f[x,x0,x1,…,xn](x−x0)(x−x1)⋯(x−xn)

为插值余项。

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

f[x,x0,x1,…,xn]=f(n+1)(ξ)(n+1)!

其中min{xi}≤ξ≤max{xi}。

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的长度不相等');endn = 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));    endend%计算牛顿插值多项式在x点上的值w = 1;y = Nt(1,1) * w;for i = 2 : n    w = w * (x - X(i-1));    y = y + Nt(i,i) * w;endend

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

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