【数值分析】牛顿插值法Matlab代码

---

前言

本文针对插值方法中的牛顿插值法进行介绍，并编写对应的Matlab代码，与先前的拉格朗日插值法类似，本文亦提供三种方法进行代码：

方法一采用多层循环进行编写，码量极小，易于复刻，但并未求出插值函数；

方法二采用符号变量结合矩阵运算，完全按照拉格朗日插值法的思路逐步编写，码量一般，但易于理解；

方法三在方法二的基础上取消了符号变量的使用，码量一般，但并未求出插值函数。

拉格朗日插值法代码详见文章：拉格朗日插值法Matlab代码

https://blog.csdn.net/Wthirteen/article/details/142589165

---

一、牛顿插值法

1.差商

在介绍牛顿插值法之前，需要对差商这一必备知识进行一个简单的学习。

对于一系列的已知点：
$$(x_i,f(x_i))i=0,1,...,n.$$
其一阶差商为：
$$f[x_i,x_{i+1}]=\frac{f(x_i)-f(x_{i+1})}{x_i-x_{i+1}}i=0,1,2,...,n-1.$$
其二阶差商为：
$$f[x_i,x_{i+1},x_{i+2}]=\frac{f[x_i,x_{i+1}]-f[x_{i+1},x_{i+2}])}{x_i-x_{i+2}}i=0,1,2,...,n-2.$$
其$m$阶差商为：
$$f[x_i,x_{i+1},....,x_{i+m}]=\frac{f[x_i,x_{i+1},....,x_{i+m-1}]-f[x_{i+1},x_{i+2},....,x_{i+m}])}{x_i-x_{i+m}}i=0,1,2,...,n-m.$$

很容易能够通过上述公式得到下表，显然，该表为一个下三角矩阵（n+1阶方阵）

函数值	一阶差商	二阶差商	…	n阶差商
$f(x_0)$
$f(x_1)$	$f[x_0,x_1]$
$f(x_2)$	$f[x_1,x_2]$	$f[x_0,x_1,x_2]$
$⋮$	$⋮$	$⋮$		$⋮$
$f(x_{n-1})$	$f[x_{n-2},x_{n-1}]$	$f[x_{n-3},x_{n-2},x_{n-1}]$
$f(x_n)$	$f[x_{n-1},x_n]$	$f[x_{n-2},x_{n-1},x_n]$		$f[x_0,x_1,...,x_n]$

2.牛顿插值函数

牛顿插值法用到了一个与拉格朗日类似的一个“算子”，但是稍有不同，故本文构造以下矩阵（n+1阶方阵）

一阶	二阶	三阶	…	$n$阶	$n+1$阶
$1$	$1$	$1$		$1$	$1$
$x-x_0$	$1$	$1$		$1$	$1$
$x-x_0$	$x-x_1$	$1$		$1$	$1$
$⋮$	$⋮$	$⋮$		$⋮$	$⋮$
$x-x_0$	$x-x_1$	$x-x_2$		$1$	$1$
$x-x_0$	$x-x_1$	$x-x_2$		$x-x_{n-1}$	$1$

将上述矩阵按行求积，即可得到一个列向量（n+1维），将其与差商表的主对角线上的元素进行线性组合，即可得到牛顿插值函数：
$$N(x)=f(x_0)+\sum _{i=1}^n\left(f[x_0,x_1,...,x_i]\cdot \prod _{j=0}^{i-1}(x-x_j)\right)$$

二、代码编写

1.方法一

代码如下：

function y=newton(X,Y,x)
%   X为已知数据点的x坐标
%   Y为已知数据点的y坐标
%   x为插值点的x坐标
%   y为各插值点函数值（以上均为列向量）
n=length(X); m=length(x);
y=zeros(m,1);
for t=1:m
    A=zeros(n,n);A(:,1)=Y;
    % 矩阵A存放的是差商表
    for  j=2:n
       for i=j:n
           % 差商表的主对角线以上部分均为0，所以i从j开始到n结束
           % 实现主对角线及以下内容的赋值
           A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1));
           % 计算差商（主对角线）
       end
    end
    for k=1:n
        p=1.0;
        for j=1:k-1
            p=p*(x(t)-X(j));
        end
        y(t)=y(t)+A(k,k)*p;        
    end
end   
end

(该段程序来源于网络，侵私删)

2.方法二

代码如下：

function y=newton(x_i,y_i,x)

%% 进行行向量检测
px=isrow(x_i);
if ~px
    x_i=x_i.';
end
py=isrow(y_i);
if ~py
    y_i=y_i.';
end
pxx=isrow(x);
if ~pxx
    x=x.';
end

%% 求解差商表
n=length(x_i);
% 初始化差商表
list_cs=zeros(n);
% 第一列为已知点函数值
list_cs(:,1)=y_i.';
list_frac2=x_i.';
% 差商
for i=2:n
    frac_1=list_cs(i-1:end-1,i-1)-list_cs(i:end,i-1);
    frac_2=-list_frac2(i:end,1)+list_frac2(1:end-i+1,1);
    list_cs(i:end,i)=frac_1./frac_2;
end

%% 求解另一部分
list_sz=ones(n,n);
for i=1:n
    list_sz(i:end,i)=x_i(i);
end

syms X
list_sz=X-list_sz;

idx=logical(triu(ones(n,n)));
list_sz(idx)=1;

list_sz=prod(list_sz,2);

%% 计算牛顿插值函数
idy1=logical(triu(ones(n,n)));
idy2=logical(triu(ones(n,n),1));
idy=xor(idy1,idy2);

Coefficient=list_cs(idy);

newton_funs=Coefficient.*list_sz;

newton_fun=sum(newton_funs);

fun=simplify(newton_fun);
disp(fun)

y=subs(fun,x);
end

该方法主要依靠于Matlab中的符号变量，计算过程中已经求解出已知数据点对应的插值函数，存储在变量fun中，但是符号变量及其表达式的运行效率相对较低。

示例函数为：
$$f(x)=\frac{1}{1+x^2}$$
选取自变量为1到10的十个整数，带入上述函数，得到因变量，以这两组数据作为插值已知数据点，求解拉格朗日插值函数，并随机生成数个自变量，利用插值函数求其因变量值，测试代码如下

clear

tic

close all

x=1:1:10;

fun=@(x)1./(1+x.^2);

y=fun(x);

plot(x,y,'r*')

x_0=sort(rand(1,20)*9);

y_0=newton(x,y,x_0);

hold on

plot(x_0,y_0,'o','Color','b')

legend("已知数据点","利用插值函数求出的点")

toc

求解结果如下图所示：

求解结果很不错。

3.方法三

代码如下：

function y=newton(x_i,y_i,x)

%% 进行行向量检测

px=isrow(x_i);
if ~px
    x_i=x_i.';
end
py=isrow(y_i);
if ~py
    y_i=y_i.';
end
pxx=isrow(x);
if ~pxx
    x=x.';
end

%% 求解差商表

n=length(x_i);
% 初始化差商表
list_cs=zeros(n);
%list_frac2=list_cs;
% 第一列为已知点函数值
list_cs(:,1)=y_i.';
list_frac2=x_i.';
% 差商
for i=2:n
    frac_1=list_cs(i-1:end-1,i-1)-list_cs(i:end,i-1);
    frac_2=-list_frac2(i:end,1)+list_frac2(1:end-i+1,1);
    list_cs(i:end,i)=frac_1./frac_2;
end

%% 求解与拉格朗日类似的算子

% 将x扩充为三维形态，便于后续操作
xx=repmat(x,length(x_i)^2,1);
xxx=reshape(xx,length(x_i),length(x_i),length(x));
first1=repmat(x_i,length(x_i),1);

% 将每一个x-x_i用矩阵表示出来，此时仍是三维矩阵
second2=xxx-first1;

idx=logical(triu(ones(n,n)));
idx2=repmat(idx,1,length(x));

idx3=reshape(idx2,n,n,length(x));

second2(idx3)=1;

list_sz=prod(second2,2);
list_sz=reshape(list_sz,n,length(x));


idy1=logical(triu(ones(n,n)));
idy2=logical(triu(ones(n,n),1));
idy=xor(idy1,idy2);
Coefficient=list_cs(idy);

%% 最终结果
newton_funs=Coefficient.*list_sz;
y=sum(newton_funs,1);

end

该方法在方法二的基础上取消了符号变量的使用，直接求解各个x对应的插值函数值y。

测试代码与方法二中测试代码相同，故仅展示求解结果：

求解结果亦很不错。

---

三、总结

方法二、三均很好的利用了matlab中矩阵的特性，从某种程度上来说，方法三与方法一其实大差不差，只有方法二中，采用了符号变量，但是一般情况来讲，符号函数的运行效率并没有想象中的高，使用的时候还需细细斟酌。方法三采用了高维矩阵，复刻代码时计算需更谨慎一些。

---

制作不易，如果这篇文章对您有所帮助，还请点个免费的赞吧，感谢！

---

参考文献

[1] 李庆扬,王能超,易大义.数值分析(第五版)[M].北京:清华大学出版社,2008:23-29.