|||
拉格朗日插值多项式,公式简洁,理论分析方便。但缺点是当改变节点或增加节点时,必须全部重新计算,这在实际计算中不方便。为此,介绍另一种形式的插值多项式 — 牛顿(Newton) 插值多项式。
一、差商及其性质
设函数 $f(x)$ 在互异点 $x_0,x_1,\cdots$ 上的值为 $f(x_0),f(x_1),\cdots$ ,定义
$f(x)$ 在 $x_i,x_j$ 上的一阶差商为
$$f[x_i,x_j] = \frac{f(x_j)-f(x_i)}{x_j-x_i},$$
$f(x)$ 在$x_i,x_j,x_k$ 上的二阶差商为
$$f[x_i,x_j,x_k] = \frac{f[x_j,x_k]-f[x_i,x_j]}{x_k-x_i},$$
递推下去,$f(x)$ 在 $x_0,x_1,\cdots,x_k$ 上的 $k$ 阶差商为
$$f[x_0,x_1,\cdots,x_k] = \frac{f[x_1,x_2,\cdots,x_k]-f[x_0,x_1,\cdots,x_{k-1}]}{x_k-x_0},$$
注:规定 $f(x)$ 在 $x_i$ 的零阶差商为$f[x_i] = f(x_i)$。
差商性质 1、$k$ 阶差商可以表示成 $k+1$ 个节点上函数值的线性组合,即
$$f[x_0,x_1,\cdots,x_k] = \sum_{j=0}^{k}\dfrac{1}{\overline{\omega}_{k+1}^\prime} f(x_j);$$
2、 差商对其节点具有对称性,即
$$f[x_{i0},x_{i1},\cdots,x_{ik}] = f[x_0,x_1,\cdots,x_k],$$
其中 $x_{i0},x_{i1},\cdots,x_{ik}$ 是~$x_0,x_1,\cdots,x_k$ 的任一排列. 这一性质由性质1直接得到.
3、若 $f(x)$ 具有 $k$ 阶连续导数,则$$f[x_0,x_1,\cdots,x_k]=\frac{f^{(k)}(\xi)}{k!},$$
其中 $\xi$ 在 $k+1$ 个节点之间.
二、牛顿插值多项式及其余项
由差商定义,把 $x$ 看成~$[a,\,b]$ 内一点,可得
$$f(x) =f(x_0)+(x-x_0)f[x_0,x],$$
由性质(2)中节点对称性, $f[x_0,x_1,x] = f[x_1,x_0,x]=\dfrac{f[x_0,x]-f[x_1,x_0]}{x-x_1}=\dfrac{f[x_0,x]-f[x_0,x_1]}{x-x_1},$
$$f[x_0,x] = f[x_0,x_1]+(x-x_1) f[x_0,x_1,x],$$
$$f[x_0,x_1,x] =f[x_0,x_1,x_2]+(x-x_2) f[x_0,x_1,x_2,x],$$
$$\cdots \cdots$$
$$f[x_0,x_1,\cdots,x_{n-1},x] =f[x_0,x_1,\cdots,x_{n-1},x_n] +(x-x_n)f[x_0,x_1,\cdots,x_{n},x],$$
只要把后一式代入前一式可得
$$f(x) = f(x_0)+(x-x_0)f[x_0,x_1]+ (x-x_0)(x-x_1)f[x_0,x_1,x_2]+\cdots+(x-x_0)(x-x_1)\cdots(x-x_{n-1})f[x_0,x_1,\cdots,x_n]+ +(x-x_0)(x-x_1)\cdots(x-x_{n})f[x_0,x_1,\cdots,x_{n},x],$$
若记
$$N_n(x) = f(x_0)+(x-x_0)f[x_0,x_1]+ (x-x_0)(x-x_1)f[x_0,x_1,x_2]+\cdots+(x-x_0)(x-x_1)\cdots(x-x_{n-1})f[x_0,x_1,\cdots,x_n],$$
$$R_n(x) =(x-x_0)(x-x_1)\cdots(x-x_{n})f[x_0,x_1,\cdots,x_{n},x],$$
则
$$f(x) = N_n(x)+R_n(x),$$
下面给出Newton基本插值的MATLAB程序。
function yi = NewtonInterp(x,y,xi)
% Newton 基本插值公式,调用格式为
% yi = NewtonInterp(x,y,xi)
% 其中
% x 为插值节点,y为节点处函数值,
% xi 为为估计函数自变量,yi 为xi处函数估计值
%
%计算差商表
n = length(x); Y= zeros(n); Y(:,1) = y';
for k = 1 : n-1
for i = 1 : n - k
if abs(x(i+k) - x(i))<eps
error('% 输入的插值节点必须互异!');
end
Y(i,k+1) = (Y(i+1,k) - Y(i,k))/(x(i+k) - x(i));
end
end
m = length(xi); yi = zeros(1,m);
for i = 1 : n % 计算Newton插值公式
z = ones(1,m);
for k = 1 : i - 1
z = z .*(xi - x(k));
end
yi = yi + Y(1,i) * z;
end
测试代码
x = [1 4 9]; y = [1 2 3]; xi = [5 6];
%x = [-2 , 0,1 , 2]; y = [17,1,2,17];
yi = NewtonInterp(x,y,xi)
>> NewtonInterpmain
yi =
2.2667 2.5000
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-7-28 00:34
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社