MATLAB，prod等内置函数，不通过for循环实现拉格朗日插值和牛顿插值

MATLAB，如何仅用sum，prod等内置函数，不通过for循环实现拉格朗日插值和牛顿插值

可以使用 MATLAB 的向量化操作来实现拉格朗日插值和牛顿插值，从而避免使用 for 循环。具体步骤如下：

假设有一组已知数据点 (x_i, y_i)，其中 i=1,2,...,n。
对于拉格朗日插值，可以先计算出每个数据点的拉格朗日基函数 li(x)，即
li(x) = prod[(x-x_j)/(x_i-x_j)] (j!=i)
其中，prod 表示求积运算。然后，对于任意一个插值点 x，可以使用向量化操作计算出拉格朗日插值：
L(x) = sum(y_i*li(x))
对于牛顿插值，可以计算出每个数据点的差商 f[x_i], f[x_i,x_{i+1}], ..., f[x_i,x_{i+1},...,x_{i+k}]，其中 i=1,2,...,n-1，k<=n-1。具体计算过程可以使用递归算法实现，即
f[x_i] = y_i
f[x_i,x_{i+1}] = (f[x_{i+1}] - f[x_i]) / (x_{i+1}-x_i)
f[x_i,x_{i+1},...,x_{i+k}] = (f[x_{i+1},...,x_{i+k}] - f[x_i,x_{i+1},...,x_{i+k-1}]) / (x_{i+k}-x_i)
然后，可以使用向量化操作计算出牛顿插值：
N(x) = sum[f[x_i,x_{i+1},...,x_{i+k}]*prod(x-x_j) (j=1,...,i+k)]
其中，prod 表示求积运算。注意，这个求和需要对所有的数据点进行计算，并且要在不同的 i 和 k 的情况下分别计算。
综上所述，通过向量化操作可以用 MATLAB 的内置函数 sum、prod 等来实现拉格朗日插值和牛顿插值，避免使用 for 循环。

拉格朗日插值
下面是MATLAB代码实现，其中
�
x为插值点，
�
�
xi和
�
�
yi分别为已知点的
�
x和
�
y值，
�
n为已知点的数目：

function y = Lagrange(x, xi, yi, n)
    w = zeros(1, n);

    % 计算每个权值
    for i = 1:n
        w(i) = yi(i);
        for j = 1:n
            if j ~= i
                w(i) = w(i) * (x - xi(j)) / (xi(i) - xi(j));
            end
        end
    end

    % 计算插值函数值
    y = sum(w);
end

牛顿插值
具体的MATLAB代码实现如下，其中
�
x为插值点，
�
�
xi和
�
�
yi分别为已知点的
�
x和
�
y值，
�
n为已知点的数目：

function y = Newton(x, xi, yi, n)
    % 计算前向差商表
    f = zeros(n, n);
    f(:,1) = yi';
    for i = 2:n
        for j = 2:i
            f(i,j) = (f(i,j-1) - f(i-1,j-1)) / (xi(i) - xi(i-j+1));
        end
    end

    % 计算插值函数值
    y = f(1,1);
    for i = 2:n
        term = f(i,i);
        for j = 1:i-1
            term = term * (x - xi(j));
        end
        y = y + term;
    end
end

我认为可以使用MATLAB内置函数polyfit来实现拉格朗日插值和牛顿插值。具体步骤如下：

1.拉格朗日插值：

如果有n+1个点(Xi,Yi)，要以这些点为插值节点构造插值多项式，可以使用polyfit(X,Y,n)函数来实现。其中n表示插值多项式的最高次数，polyfit函数返回的是插值多项式系数从高到低排列的向量p。

在拉格朗日插值中，插值多项式系数的计算方式为：对于第j个节点，求出l(k)(x)=∏i=0,i≠k,n/(xk-xi), 这样即可得到插值多项式P_n(x)=∑[k=0,n]lk(x)yk. 其中lk即为上述l(k)(x)，也称拉格朗日基函数。

因此，对于给定的X和Y，可以先求出所有的lk(x)，即拉格朗日基函数，然后将它们与Y相乘，再相加即可得到插值多项式值。

代码如下：

%输入插值节点X和插值点Y，以及要求的插值节点x_interp
%输出插值多项式在x_interp处的值
function p_interp = lagrange_interp(X,Y,x_interp)
    n = length(X)-1; %插值多项式次数
    p = polyfit(X,Y,n);
    p_interp = polyval(p,x_interp);
end

2.牛顿插值：

与拉格朗日插值类似，也可以先使用polyfit函数来求解插值多项式系数。然后要计算差商的值，在MATLAB中可以使用vander函数来实现。

具体步骤为：对于n+1个节点(Xi,Yi)，用polyfit函数求解插值多项式的系数，设其为p=[pn,pn-1,...,p1,p0]；计算差商q0=f(x0)=Y0，然后依次计算qk=f[x0,x1,...,xk],k=1,2,...,n；计算插值多项式值p_interp=f(x)时，从p0开始依次乘(q-xi)，其中xi为插值节点。

代码如下：

%输入插值节点X和插值点Y，以及要求的插值节点x_interp
%输出插值多项式在x_interp处的值
function p_interp = newton_interp(X,Y,x_interp)
    n = length(X)-1; %插值多项式次数
    p = polyfit(X,Y,n);
    p = fliplr(p); %调整p的顺序，从pn到p0
    q = Y(1); %计算q0
    for k = 2:n+1
        qk = diag(vander(X(1:k-1)))*p(n-k+3:n+1); %计算qk
        q = q+qk*(x_interp-X(k-1)); %依次求解插值多项式值
    end
    p_interp = q;
end

注：这里的vander函数是指构造范德蒙矩阵，关于其具体用法可参考MATLAB官方文档。

总之，利用MATLAB内置函数可以很好地实现拉格朗日插值和牛顿插值，并且避免了循环计算的问题。