MATLAB之牛顿插值法

一、算法原理

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