一、算法原理
1、问题引入
上一篇博客我们介绍了拉格朗日插值法,我们现在来回顾其定义:
对于插值区间【a,b】上一系列插值节点x0,x1,x2,......xn及其函数值y0、y1、y2......yn;构造一个简单的简单易算的近似函数g(x)≈f(x),满足条件g(xi)=f(xi) (i=0,1,2...,n)。
其拉格朗日插值多项式的形式为:
g(x)=L0(x0)*y0+L1(x1)*y1+......+Ln(xn)*yn,
其中Li(xi)=[(x-x0)*(x-x1)*......*(x-x(i-1))*(x-x(i+1))*.....*(x-xn)] / [(xi-x0)*(xi-x1)*......*(xi-x(i-1))*(xi-x(i+1))*.....*(xi-xn)]。
优点:具有严格的规律性,便于记忆。
缺点:计算量大、不具有承袭性。
2、牛顿插值法
Lagrange 插值虽然易算,但若要增加一个节点时, 全部基函数Li(x) 都需要重新计算。能否改写其形式,使其每增加一个点,只附加一项(重新计算一项)。
设Nn(x)=A0+A1*(X-X0)+A2*(X-X0)(X-X1)+......+An(X-X0)*(X-X1)......(X-X(n-1))
利用Nn(xi)=f(xi),写成矩阵形式:
1 0 ...... 0 [A0 [f(x0) 1 x1-x0 ...... 0 A1 f(x1) 1 x2-x0 (x2-x0) (x2-x1) ...... 0 * ...... = ...... ... ... ...... ...... ... ... ... ...... ...... 1 xn-x0 (xn-x0) (xn-x1) ...... (xn-x0)*...(xn-x(n-1)) An] f(xn)] 根据该线形方程组可以求得:A0=f(x0),A1=[f(x1)-f(x0)]/[x1-x0], A2={[f(x2)-f(x0)]/[x2-x0]-f(x1)-f(x0)]/[x1-x0]} / (x2-x1)
|
3差商的定义
设已知函数f (x)在一系列互不相等的节点x0,x1,x2,......xn,(即在i≠j时,xi≠xj)上的函数值 f(xi) ,
f (x)在点xi , xj处的一阶差商为:
f (x)在点xi,xj,xk处的二阶差商为:
f(x)在点x0,x1,x2,......xk的k阶差商为:
规律总结:高阶差商是两个低一阶差商的差商。
利用差商的定义,可得Nn(x)的系数Aj :
从而可得牛顿插值多项式:
这样就实现了,每增加一个节点,多项式只增加一项,克服了拉格朗日插值多项式的缺点。
差商可以通过差商表计算。
4、牛顿插值多项式
了解了差商,我们就可以利用差商来表示牛顿插值多项式:
二、matlab程序
%% 牛顿插值 (创建函数脚本方可运行)
close
clc
clear
x=[2 3 5 6 9];
y=[4 8 12 15 23];
m=length(X);
c=cha(X,Y);
syms x
for i=1:m-1
D1=x-X(1:i);
D(i)=prod(D1);
end
N=sym2poly(expand(sum(D.*c)+Y(1)))
x=X(1):0.1:X(end);
y=polyval(N,x);
plot(x,y)
hold on
plot(X,Y,'o')
function C=cha(X,Y)
% C存储差商结果
% X是横坐标
m=length(X);
d=diff(Y)./diff(X); %一阶差商
C(1)=d(1); %取出第一个一阶差商做系数
for i=2:m-1 %i阶差商 依次算出C(2)、C(3)...C(m-1)阶茶商 ,储存在数组C中
X1=X(i+1:end);
n=length(X1);
X2=X(1:n);
d=diff(d)./(X1-X2);
C(i)=d(1);
end
end