Matlab主程序

X=[0.2,0.4,0.6,0.8,1.0];
Y=[0.98,0.92,0.81,0.64,0.38];
M=1;
m=[0,1,11,10];
a=length(m);
x=zeros(size(m));%生成和m相同的零矩阵
for i=1:a
    x(i)=0.2+0.08*m(i);
end
x1=linspace(0.2,1.0,10);
[y R A]=Newton(X,Y,x1,M);
[y1 R1 A1]=Newton(X,Y,x,M);
s=interp1(X,Y,x1,'spline');
s1=interp1(X,Y,x,'spline');
plot(X,Y,'or',x1,y,'-b',x1,s,'-k',x,y1,'sb',x,s1,'>k');
legend('Raw','Newton line','Spline line','Newton dots','Spline dots');

Newton插值法函数

function [y,R,A]= Newton(X,Y,x,M)%其中A是如书中31页所示的均差表，理清楚其中的数字关系，就能完成编程
n = length(X); m = length(x);
for t=1:m
    z=x(t);
    A=zeros(n,n);
    A(:,1)=Y';
    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
    s=A(1,1);
    for k=2:n
        C=poly(X(1:k-1));
        s=s+A(k,k)*polyval(C,z);%虽然看上去思路清晰了，但是感觉不是效率最高的，没有用到比较高级的数学工具。
    end
    y(t)=s;
    R(t)=M*q1/c1;
end

结果图

1.jpg

重点说明

poly(a1,a2,...,an)表示的是(x-a1)(x-a2)...(x-an)的多项式，从高次到常数项的系数矩阵。
polyval(A,z)表示的是，以A行列式为系数的多项式，x=z时的取值。
plot函数中，单引号中的语句可以调整最后呈现的方式。可以通过在主窗口输入help plot进行查阅。
interp1(X,Y,x,'method')是matlab自带的插值函数。其中X,Y表示插值点，x表示自变量，‘method’共有一下几种插值方式：

      'nearest'  - nearest neighbor interpolation
      'linear'   - linear interpolation
      'spline'   - piecewise cubic spline interpolation (SPLINE)
      'pchip'    - shape-preserving piecewise cubic interpolation
      'cubic'    - same as 'pchip'
      'v5cubic'  - the cubic interpolation from MATLAB 5, which does not
                   extrapolate and uses 'spline' if X is not equally
                   spaced.