牛顿多项式插值函数实现问询及解法对比请求

原需求说明 #

编写函数，输入自变量x向量、对应因变量y向量、待插值向量，输出牛顿多项式系数向量及插值结果，代码需添加详细注释。

我的尝试代码 #

function yint = Newtint(x,y,xx)
    n = length(x);

    if length(y)~=n, error('x and y must be same length'); end
    
    b = zeros(n,n);
    b(:,1) = y(:); % the (:) ensures that y is a column vector.
    
    for j = 2:n
        for i = 1:n-j+1
            b(i,j) = (b(i+1,j-1)-b(i,j-1))/(x(i+j-1)-x(i));
        end
    end
    
    xt = 1;
    yint = b(1,1);
    
    for j = 1:n-1
        xt = xt*(xx-x(j));
        yint = yint+b(1,j+1)*xt;
    end
    
    
    % input:
    % x = independent variable
    % y = dependent variable
    % xx = value of independent variable at which
    % interpolation is calculated
    % output:
    % yint = interpolated value of dependent variable
    % compute the finite divided differences in the form of a
    % difference table

尝试代码的不足 #

• 仅输出插值结果yint，未按要求输出牛顿多项式的系数向量
• 未处理待插值向量xx为多维的情况，当xx是向量时计算逻辑会出错
• 注释位置零散，未嵌入代码流程，可读性不足

正确实现代码 #

function [coeff, yint] = NewtonInterp(x, y, xx)
    % 牛顿多项式插值函数
    % 输入参数:
    %   x       - 自变量样本向量（长度为n，对应n个插值节点）
    %   y       - 因变量样本向量（需与x长度一致）
    %   xx      - 待插值的自变量向量（支持单个值或多个值）
    % 输出参数:
    %   coeff   - 牛顿多项式系数向量，形式为[a1, a2, ..., an]，对应多项式：
    %             P(x) = a1 + a2*(x-x1) + a3*(x-x1)(x-x2) + ... + an*(x-x1)...(x-x(n-1))
    %   yint    - 对应xx的插值结果向量，与xx长度一致

    % 获取节点数量
    n = length(x);
    % 检查x与y长度是否匹配
    if length(y) ~= n
        error('自变量x和因变量y的长度必须一致！');
    end

    % 初始化差商表，n行n列，第一列存储y值
    diff_table = zeros(n, n);
    diff_table(:, 1) = y(:);  % 将y转为列向量存入第一列

    % 计算各阶差商
    for j = 2:n  % j表示当前计算的差商阶数（从2阶到n阶）
        for i = 1:n - j + 1  % 计算第j阶差商的第i个值
            % 差商公式：[xi,...,xi+j-1] = ([xi+1,...,xi+j-1] - [xi,...,xi+j-2]) / (xi+j-1 - xi)
            diff_table(i, j) = (diff_table(i+1, j-1) - diff_table(i, j-1)) / (x(i+j-1) - x(i));
        end
    end

    % 提取牛顿多项式系数：差商表的第一行
    coeff = diff_table(1, :);

    % 计算插值结果，支持xx为向量的情况
    m = length(xx);
    yint = zeros(m, 1);  % 初始化插值结果向量

    for k = 1:m  % 遍历每个待插值点
        x_current = xx(k);
        % 初始化插值多项式值与累积项
        p_val = coeff(1);
        term = 1;
        for j = 2:n  % 逐步累加各阶项
            term = term * (x_current - x(j-1));  % 累积(x-x1)(x-x2)...(x-x(j-1))
            p_val = p_val + coeff(j) * term;
        end
        yint(k) = p_val;
    end
end

代码说明 #

• 差商表的第一行即为牛顿插值的系数向量，对应多项式的各项系数
• 新增了对向量型待插值输入xx的支持，每个插值点独立计算
• 注释嵌入代码流程中，清晰说明每一步的作用
• 严格区分输出参数，同时返回系数向量和插值结果

内容的提问来源于stack exchange，提问作者Daisy Ridley