多项式插值
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 直接法
根据定理一,假设插值多项式为
由插值条件 yi=P(xi) ( i=0,1,…,n ),我们得到关于系数 an,an−1,…,a1,a0 的线性方程组
通过求解这个线性方程组,即得到插值多项式。
优点:直接,性质一目了然。
缺点:待求解的线性方程组的系数矩阵为 范德蒙德 (Vandermonde)矩阵,它是一个病态矩阵,这使得在实际求解方程组时将产生很大的误差。
3.2 拉格朗日多项式插值
拉格朗日(Lagrange)多项式插值的原理是:先构造一组 拉格朗日基函数
Li(x)
(
i=0,1,…,n
),这些基函数为次数不超过
n
的多项式,且具有性质
然后将这些基函数做线性组合,得到拉格朗日插值多项式
容易验证,多项式 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的长度不相等');
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) 关于点 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 阶 | 2 阶 | ⋯ | 因子 | ||
---|---|---|---|---|---|---|
x0 | f(x0) | 1 | ||||
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,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的长度不相等');
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 拉格朗日多项式插值和牛顿多项式插值的比较
拉格朗日多项式插值的计算量大于牛顿多项式插值的计算量。特别地,当新增一个插值点时,拉格朗日插值需要重新计算全部的基函数,而牛顿插值只需计算均差表中新的一行的值即可。