newton.m

《MATLAB实用指南》系列丛书DE源代码

function s=Newton(x,y,x0,nn)
  % Newton插值，x与y为已知的插值点及其函数值
  % x0为要求的插值点的x坐标值。nn为Newton插值多项式的阶数
  nx=length(x);
  ny=length(y);
  if nx~=ny
      warning('向量x与y的长度应该相同')
      return
  end
  m=length(x0);
  %按照公式，对要求的插值点x0的每个元素进行计算
  for i=1:m
      t=0.0;
      j=1;
      yy=y;
      kk=j;
      %求各级均差
      while(kk<=nn)
          kk=kk+1;
          for k=kk:nx
              yy(k)=(yy(k)-yy(kk-1))/(x(k)-x(kk-1));
          end
      end
      % 求差值结果
      t=yy(1);
      for k=2:nn
          u=1.0;
          jj=1;
          while(jj<k)
              u=u*(x0(i)-x(jj));
              jj=jj+1;
          end
          t=t+yy(k)*u;
      end
      s(i)=t;
  end