1）Newton插值公式源程序：

function newtoncz(a,b,n,f)

ln=length(n);                    %给定n的个数

for k=1:ln

    m=n(k)-1;                     %等分份数

    y=zeros(1,m+1);

    A=zeros(m+1);

    w=zeros(1,m+1);

    h=(b-a)/m;  

    for i=1:m+1

        x(i)=a+(i-1)*h;

        y(i)=subs(f,findsym(f),x(i));    %插值节点函数值

    end

    A(:,1)=y';

    for i=2:m+1

        for j=i:m+1

            A(j,i)=(A(j,i-1)-A(j-1,i-1))/(x(j)-x(j-i+1));

        end

    end

    p=A(1,1);

w=vpa(w,4);

syms X;

    w(1)=X-x(1);

    for i=2:m+1

    w(i)=w(i-1)*(X-x(i));  

    p=A(i,i)*w(i-1)+p;

    p=simplify(p);

    end

   p=vpa(p,4);

   fprintf('n=%d的newton插值多项式为：',n(k));

    disp(p);

    h=ezplot(p,[a,b]);

    if k==1

        set(h,'color','r');

    elseif k==2

        set(h,'color','b');

    elseif k==3

        set(h,'color','g');

    else

        set(h,'color','m');

    end

    hold on

end

h=ezplot(f,[a,b]);

set(h,'linewidth,2,'color,'k');

2)在命令窗口输入a，b,n,f,得到Newton插值公式以及原函数、差值函数的图像：

>> a=-1;

>> b=1;

>>n=[5,7,10,13];

>> syms X;

>> f=1/(1+25*X.^2);

>> newtoncz(a,b,n,f)

运行结果：

n=5的newton插值多项式为：3.316*X^4-4.277*X^2+1.

n=7的newton插值多项式为：

-13.13*X^6+20.96*X^4+.3475e-14*X^3-8.784*X^2-.2420e-15*X+1.

n=10的newton插值多项式为：

.3730e-13*X^9+21.62*X^8-.6948e-13*X^7-44.92*X^6+.4163e-13*X^5+30.73*X^4+.7857e-14*X^3-8.261*X^2+.1321e-14*X+.8615

n=13的newton插值多项式为：

909.9*X^12-.3411e-12*X^11-2336.*X^10+.6632e-12*X^9+2202.*X^8-.1573e-12*X^7-955.4*X^6+.2341e-12*X^5+198.7*X^4-.3535e-13*X^3-19.58*X^2+.7881e-15*X+1.000

 ------红色（n=5）------蓝色（n=7）-----绿色（n=10）-----紫红色（n=13）-----黑色（原函数）

分析：由程序可知，当插值节点个数变化时，Newton插值多项式的结构不改变，插值多项式易于构造。只要给定插值函数，插值区间，插值节点，就可以很容易的得出插值多项式及其图像。由图像可知，随着n的增大，插值精度逐渐提高，但并不是插值节点越多逼近精度越高，因为随着节点数目的增加导致插值多项式的次数增高，而高次多项式的振荡次数增多可能使插值多项式在非节点处的误差变得很大。如当n=10和13时，图像两边与原函数偏离较大。下载本文