MATLAB牛顿插值

X为各个插值节点，Y为各阶均差

程序如下：

function [A,C,L,wcgs,Cw]= newploy(X,Y)

n=length(X); A=zeros(n,n); A(:,1)=Y';

        s=0.0; p=1.0; q=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        

             b=poly(X(j-1));q1=conv(q,b); c1=c1*j;  q=q1;

     end

 C=A(n,n); b=poly(X(n)); q1=conv(q1,b);     

for k=(n-1):-1:1

C=conv(C,poly(X(k))); d=length(C); C(d)=C(d)+A(k,k);

end

L(k,:)=poly2sym(C); Q=poly2sym(q1);

syms M

wcgs=M*Q/c1; Cw=q1/c1;

运行：书P32 例4

X=[0.40 0.55 0.65 0.80 0.90 1.05];

Y=[0.41075 0.57815 0.69675 0.88811 1.02652 1.25382];

[A,C,L,wcgs,Cw]= newploy(X,Y)

输出结果：

A =

  Columns 1 through 4

    0.4108         0         0         0

    0.5782    1.1160         0         0

    0.6967    1.1860    0.2800         0

    0.8881    1.2757    0.35    0.1973

    1.0265    1.3841    0.4335    0.2130

    1.2538    1.5153    0.5249    0.2287

  Columns 5 through 6

         0         0

         0         0

         0         0

         0         0

    0.0312         0

    0.0314    0.0003

C =

  Columns 1 through 4

    0.0003    0.0303    0.1236    0.0296

  Columns 5 through 6

    0.9901    0.0013

L =

27028192032443/9223372036854775808*x^5+2181259916290823/72057594037927936*x^4

+4453713600685277/360287970163968*x^3+2134103104007001/72057594037927936*x^2

+18190441637495/9007199254740992*x+229755220417/180********481984

wcgs =

1/720*M*(x^6-87/20*x^5+3097/400*x^4-57681/8000*x^3+4166716301494349/11259906842624*x^2

-44711044588153/4503599627370496*x+35001188126157/360287970163968)

Cw =

  Columns 1 through 4

    0.0014   -0.0060    0.0108   -0.0100

  Columns 5 through 7

0.0051   -0.0014    0.0002