牛顿插值法的 Matlab 实现

意外的发现了卷积可以用来求多项式相乘的系数。

特点

在学插值之前就已经在 Matlab 上用过多项式插值了，也就是 polyfit 和 polyval。所以在实现牛顿插值的时候，我就希望我的实现在用法上要和 polyfit 的用法是相似的。

函数 polyfit 的作用是返回多项式的系数，而函数 polyval 是根据多项式的系数，得到一系列点对应的多项式值。那么我的 Interpolation_Newton 也应该返回多项式的系数。

然后我就遇到了问题，我需要把含有商差的多项式转化为标准的多项式，这就需要做多项式乘法，但我不知道如何高效地实现（其实是懒得写这个程序）。查了别人关于牛顿插值的实现代码，发现有不少实现都用到了卷积计算(conv)。最后我在 MathWorks 的文档页面得到了我想要的东西：

w = conv(u,v) 返回向量 u 和 v 的卷积。如果 u 和 v 是多项式系数的向量，对其卷积与将这两个多项式相乘等效。

之前也看过一些关于卷积的科普，但好像都没提到卷积还和多项式相乘有关系。

代码

最后写出来就是这样：

%% INIT
% 《数值分析方法与应用》张宏伟、孟兆良编著；大连理工大学出版社
% 第 217 页：四、插值与逼近 第 5 题
clear;clc;
format short;
%% 初始化参数
% 区间
range = [0 1];
% 函数定义
func = @(x) sin(pi .* x);
%% 计算与输出
% 已知数据
xx = range(1):0.1:range(2);
yy = func(xx);
% 插值数据
p = Interpolation_Newton(xx, yy);
xx1 = range(1):0.05:range(2);
yy1 = polyval(p,xx1);
% 原函数
x = range(1):0.01:range(2);
plot(x,func(x));   hold on;
% 插值点图像
plot(xx1,yy1,'--o'); hold on;
legend('原函数', '插值')

%% 函数定义
function factors =  Interpolation_Newton(x, y)
    n = length(x);
    dd = zeros(n,n);
    % 计算均差
    dd(1,:) = y;
    for k = 1:n-1
        for r = 1:n-k
            dd(k+1, r) = ( dd(k, r+1) - dd(k, r) ) / (x(r+k)-x(r));
        end
    end
    % 计算各阶多项式系数
    p = zeros(n,n);
    p(1,end) = 1;
    for k = 2:n
        % conv (卷积计算):
        % 对两个多项式系数进行卷积计算可得到多项式相乘的系数
        p(k,:) = conv(p(k-1,2:end), [1, -x(k-1)]);
    end
    factors = sum(dd(:,1) .* p);
end