6.3.5 牛顿插值法的MATLAB综合程序
求牛顿插值多项式、差商、插值及其误差估计的MATLAB主程序
function [y,R,A,C,L]=newdscg(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(k)= polyval(C, z); end
R=M*q1/c1;L(k,:)=poly2sym(C);
例6.3.6 给出节点数据f(?4.00)?27.00，f(0.00)?1.00，f(1.00)?2.00，
f(2.00)?17.00，作三阶牛顿插值多项式，计算f(?2.345)，并估计其误差.
解 首先将名为newdscg.m的程序保存为M文件，然后在MATLAB工作窗口输入程序
>> syms M,X=[-4,0,1,2]; Y =[27,1,2,17]; x=-2.345; [y,R,A,C,P]=newdscg(X,Y,x,M)
运行后输出插值y?f(?2.345)及其误差限公式R，三阶牛顿插值多项式P及其系数向量C，差商的矩阵A如下
y =
22.3211 R =
1323077530165133/562949953421312*M（即R =2.3503*M） A=
27.0000 0 0 0 1.0000 -6.5000 0 0 2.0000 1.0000 1.5000 0 17.0000 15.0000 7.0000 0.9167 C =
0.9167 4.2500 -4.1667 1.0000 P =
11/12*x^3+17/4*x^2-25/6*x+1