牛顿插值法的 Matlab 实现

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

特点

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

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

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

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

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

代码

最后写出来就是这样:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
%% 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

作者

uint128.com

发布于

2023-01-08

更新于

2023-01-08

许可协议

评论