newpoly.m

来自「数值计算实践是我们必备的一门学科！这是我在学习之余利用MATLAB编写的程序」· M 代码 · 共 23 行

function [A,C,L,wcgs,Cw]=newpoly(X,Y)
%A为差商的矩阵;
%L为牛顿插值多项式;
%C为牛顿插值多项式的系数向量;
%wcgs为插值余项公式;
%Cw为余项和插值函数比值的系数向量.
%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;