拉格朗日函数

 function f=language(x,y,x0)
%求已知数据点的拉格朗日插值多项式
%已知数据点的横坐标：x
%已知数据点的纵坐标：y
%求得的拉格朗日多项式或在x0处的插值：f
syms t;
if (length(x)==length(y))
    n=length(x);
else
    disp('x和y的维数不同');
    return;
end
f=0.0;
for (i=1:n)
    l=y(i);
    for j=1:(i-1)
        l=l*(t-x(j))/(x(i)-x(j));
    end;
    for j=(i+1):n
        l=l*(t-x(j))/(x(i)-x(j));%计算拉格朗日基函数
    end;
    f=f+l;                       %计算拉格朗日插值函数
    simplify(f);                 %化简
    if i==n
        if nargin==3
            f=subs(f,'t',x0);    %计算插值点的函数值
        else
            f=collect(f);        %将多项式展开
            f=vpa(f,4);          %将插值多项式的系数化成4位精度的小数
        end
    end
end