利用均差的牛顿插值法(Newton)

函数f的零阶均差定义为 ，一阶定义均差为：

一般地，函数f 的k阶均差定义为：

或者上面这个式子求的k+1阶均差

利用均差的牛顿插值法多项式为：

简单计算的时候可以观看下面的差商（均差）表：

怎么利用差商表计算，可以看下面这个例子：

正常的话还有一个余项，在本文中先不考虑了。

总和上面的计算方法可以归纳出算法的大致思想：先计算差商表，类似于乘法口诀的思路，两个for循环就可以计算出，然后对于每一次内for循环以后，计算出了第一列，接着把相对应的f(x)计算出来，接着进入第二列的计算，接着计算相应的f(x).......一直到计算完毕最后一个f(x)，把所有的f(x)相加，便是最终的插值。

为了便于写算法，以五个样本点为例，计算了一下差商表：

【注】这里的覆盖是就是把这一列计算的值覆盖到前一列的对应地方，这样便于找到规律，可以发现分子始终是y(i+1)-y(i)，而分母则与列号有关系，详看代码，更易于理解

Newton1.m

function f = Newton1(x,y,x0)%求已知数据点的均差形式牛顿插值多项式%已知数据点的x坐标向量：x%已知数据点的y坐标向量：y%插值点的x坐标：x0%求得的均差形式牛顿插值多项式或x0处的插值：fsyms t;%计算输入的x y是否长度相等if(length(x)==length(y))    n=length(x);    c(1:n)=0.0;else    dis('x和y的维数不相等！');    return;  endf = y(1);%第0列的f(x)就是y(1)本身y1 = 0; %这个y1不是y(1)，存的是差商表后面的值l  = 1; %l是用来算f(x)后面对应的(x-x1)(x-x2).....的 for(i=1:n-1)       for j=1:i        y1(j)=0;    end    for(j=i+1:n)        y1(j) = (y(j)-y(j-1))/(x(j)-x(j-i)); %利用前面的列计算后面列的值    end    c(i) = y1(i+1);         l = l*(t-x(i));      f = f + c(i)*l;    simplify(f);    y = y1;        if(i==n-1)        if(nargin == 3)            f = subs(f,'t',x0);%替换函数，用后面的替换前面的，把t替换为x0        else            f = collect(f);                %将插值多项式展开            f = vpa(f, 6);        end    endend

Newton1Insert.m

% x=[1 1.2 1.8 2.5 4];% y=[1 1.44 3.24 6.25 16];% f=Newton1(x,y)% f=Newton1(x,y,2.0)% x=[0.40 0.55 0.65 0.80 0.90];% y=[0.41075 0.57815 0.69675 0.88811 1.02652];% f=Newton1(x,y,0.596)x1=0:2*pi;y1=sin(x1);xx=0:0.2:2*pi;yy=Newton1(x1,y1,xx);plot(x1,y1,'o:',xx,yy,'r')