牛顿（Newton)插值及其MATLAB程序

拉格朗日插值的优点是格式整齐和规范，有误差估计方式，它的缺点是没有承袭性，当需要增加节点时，必须重新计算插值的基函数li(x).本文给出具有承袭性的牛顿插值法及其MATLAB程序。与牛顿插值有关的差商相关概念及其性质请参看《计算方法－差商的概念及性质》

一、牛顿插值公式



二、牛顿插值法的MATLAB的综合程序



编写M文件

%newton.m
%求牛顿插值多项式、差商、插值及其误差估计的MATLAB主程序
%输入的量:X是n+1个节点(x_i,y_i)(i = 1,2, ... , n+1)横坐标向量，Y是纵坐标向量，
%x是以向量形式输入的m个插值点，M在[a,b]上满足｜f~(n+1)(x)｜≤M
%注：f~(n+1)(x)表示f(x)的n+1阶导数
%输出的量：向量y是向量x处的插值，误差限R，n次牛顿插值多项式L及其系数向量C，
%差商的矩阵A
function[y,R,A,C,L] = newton(X,Y,x,M)
n = length(X);
m = length(x);
for t = 1 : m
z = x(t);
A = zeros(n,n);
A(:,1) = Y';
s = 0.0; p = 1.0; q1 = 1.0; c1 = 1.0;
for j = 2 : n
for i = j : n
A(i,j) = (A(i,j-1) - A(i-1,j-1))/(X(i)-X(i-j+1));
end
q1 = abs(q1*(z-X(j-1)));
c1 = c1 * j;
end
C = A(n, n); q1 = abs(q1*(z-X(n)));
for k = (n-1):-1:1
C = conv(C, poly(X(k)));
d = length(C);
C(d) = C(d) + A(k,k);%在最后一维，也就是常数项加上新的差商
end
y(t) = polyval(C,z);
R(t) = M * q1 / c1;
end
L = poly2sym(C);


例：用Lagrange插值来求sinx在某点的值，并估计其误差，已知sin0°  = 0, sin30°  = 0.5, sin45° = 0.7071, sin60° = 0.8660, sin90° = 1.

X                       0                         pi/6                         pi/4                        pi/3                     pi/2

Y                       0                         0.5                        0.7071                  0.8660                   1 

解：

>> X = [0 pi/6 pi/4 pi/3 pi/2];
>> Y = [0 0.5 0.7071 0.8660 1];
>> x = linspace(0,pi,50);
>> M = 1;
>> [y,R,A,C,L] = newton(X, Y, x, M);
>> y1 = sin(x);
>> errorbar(x,y,R,'.g')
>> hold on
>> plot(X, Y, 'or', x, y, '.k', x, y1, '-b');
>> legend('误差','样本点','牛顿插值估算','sin(x)');