函数逼近——(Newton)牛顿插值法 | 北太天元

Lagrange 插值

• 优点：容易得到，公式紧凑

• 不足：当增删节点时，计算需重新进行

解决方案：采用新的基函数，逐次生成多项式。 下面介绍 Newton 插值

一、 Newton插值

已知 $n+1$ 个不同的样本点 $(x_i,y_i),i=0,1,2,\dots ,n$
$n$ 次Newton插值多项式
$$N_n(x)=C_0{\varphi }_0(x)+C_1{\varphi }_1(x)+\cdots +C_n{\varphi }_n(x)$$
其使用的基函数如下：

• ${\varphi }_0(x)=1$
• ${\varphi }_1(x)=(x-x_0)$
• ${\varphi }_2(x)=(x-x_0)(x-x_1)$
• $⋮$
• ${\varphi }_n(x)=(x-x_0)(x-x_1)\cdots (x-x_{n-1})$

基函数应满足

$$\left(\begin{array}{cccc}{\varphi }_0(x_0)& 0& \cdots & 0\\ {\varphi }_0(x_1)& {\varphi }_1(x_1)& \cdots & 0\\ ⋮& ⋮& & ⋮\\ {\varphi }_0(x_n)& {\varphi }_1(x_n)& \cdots & {\varphi }_n(x_n)\end{array}\right)\left(\begin{array}{c}{C}_0\\ C_1\\ ⋮\\ C_n\end{array}\right)=\left(\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y_n\end{array}\right)$$

解这个方程组，可以得到 $C_0,C_1,\dots ,C_n$

结合差商的定义

$$C_0=f(x_0),C_1=f[x_0,x_1],C_2=f[x_0,x_1,x_2],\dots ,C_n=f[x_0,x_1,x_2,\dots ,x_n]$$

得

$$N_n(x)=f(x_0)+f[x_0,x_1](x-x_0)+\cdots +f[x_0,x_1,\cdots ,x_n](x-x_0)(x-x_1)\cdots (x-x_{n-1})$$

其中，差商可由差商表进行计算,

$$\begin{array}{cccccc}{x}_i& f[x_i]& f[x_i,x_{i+1}]& f[x_i,x_{i+1},x_{i+2}]& f[x_i,x_{i+1},x_{i+2},x_{i+3}]& \cdots \\ x_0& f[x_0]& & & & \\ x_1& f[x_1]& f[x_0,x_1]& & & \\ x_2& f[x_2]& f[x_1,x_2]& f[x_0,x_1,x_2]& & \\ x_3& f[x_3]& f[x_2,x_3]& f[x_1,x_2,x_3]& f[x_0,x_1,x_2,x_3]& \\ \cdots & & & & & \end{array}$$

二、算法的简单介绍

💖 Newton 插值 (这里只强调差商表对计算的节省)

输入:

• $n+1$ 个样本点的横纵坐标分别构成的向量 $x0,y0$
• 目标近似点 $x$

输出:

• 近似点的值 $y$

实现步骤

1. 步骤 $1$: 计算差商表,用矩阵$D$表示,有 $n$ 行,$n+1$ 列

$$\begin{array}{ccccc}{x}_0& f[x_0]& & & \\ x_1& f[x_1]& f[x_0,x_1]& & \\ x_2& f[x_2]& f[x_1,x_2]& f[x_0,x_1,x_2]& \\ x_3& f[x_3]& f[x_2,x_3]& f[x_1,x_2,x_3]& f[x_0,x_1,x_2,x_3]\\ \cdots & & & & \end{array}$$

2. 步骤 $2$: 计算基函数 ${\varphi }_i(x)$
3. 步骤 $3$:
   $$N_n(x)=f(x_0)+f[x_0,x_1](x-x_0)+\cdots +f[x_0,x_1,\cdots ,x_n](x-x_0)(x-x_1)\cdots (x-x_{n-1})$$

三、北太天元源程序

计算差商表

function [D] = divided_differences(x0,y0,c)
% Newton插值中用到的差商计算
% Input: 互异样本点构成的向量 (横) x0, (纵) y0， c:控制列循环的次数
% Output: 差商表 D
%   Version:            1.0
%   last modified:      04/01/2024
n = length(x0);
D = zeros(n,n+1); 
D(:,1) = x0 ;D(:,2) = y0;
if nargin < 3
    c = n+1;
end
% 列循环
for k = 3:1:c  
    % 行循环
    for i =k-1:1:n 
        D(i,k) = D(i,k-1) - D(i-1,k-1);
        D(i,k) = D(i,k)/(x0(i)-x0(i-(k-2)));
        %D   % 取消注释 D 可以观察每次的变化
    end
end
end

将上述代码保存为 divided_differences.m 文件。

Newton插值

function [N,D] = Newton_interp(x0,y0,x)
% Newton 插值法 插值逼近
% Input: 样本点构成的向量x0 (升序), (纵) y0 
%        x0,y0 是行向量
%        x  
% Output: N 插值结果 , D 差商表
%   Version:            1.0
%   last modified:      04/01/2024
n = length(x0);
D = divided_differences(x0,y0);
D_need = diag(D(:,2:n+1)); % 列
% 表示基函数
n1 = length(x);
for k = 1:1:n1
    xjx = [1,x(k)-x0(1:n-1)]; 
    for i = 2:length(xjx)
        xjx(i) = xjx(i)*xjx(i-1);
    end
    N(k) = xjx * D_need;
end
end

将上述代码保存为 Newton_interp.m 文件。

四、简单实现一下

例子1

%% 测试差商计算 
clc;clear;
x0 = [0.4 0.5 0.6 0.7 0.8];
y0 = [-0.916291 -0.693147 -0.510826 -0.357765 -0.223143];
D =divided_differences(x0,y0,4)

其中计算差商表函数中的 c ，控制着算到第 c-2 阶差商

所得差商表，第1列是 x0，第2列 为 y0，第3列为一阶差商，第4列为二阶差商

c = 3 时

D =

0.4000   -0.9163    0.0000    0.0000    0.0000    0.0000
0.5000   -0.6931    2.2314    0.0000    0.0000    0.0000
0.6000   -0.5108    1.8232    0.0000    0.0000    0.0000
0.7000   -0.3578    1.5306    0.0000    0.0000    0.0000
0.8000   -0.2231    1.3462    0.0000    0.0000    0.0000

c = 4时：

D =

0.4000   -0.9163    0.0000    0.0000    0.0000    0.0000
0.5000   -0.6931    2.2314    0.0000    0.0000    0.0000
0.6000   -0.5108    1.8232   -2.0412    0.0000    0.0000
0.7000   -0.3578    1.5306   -1.4630    0.0000    0.0000
0.8000   -0.2231    1.3462   -0.9220    0.0000    0.0000

c 不填时

D =

0.4000   -0.9163    0.0000    0.0000    0.0000    0.0000
0.5000   -0.6931    2.2314    0.0000    0.0000    0.0000
0.6000   -0.5108    1.8232   -2.0412    0.0000    0.0000
0.7000   -0.3578    1.5306   -1.4630    1.9272    0.0000
0.8000   -0.2231    1.3462   -0.9220    1.8035   -0.3092

通过c的控制，我们可以只算到我们想要的第几阶差商。

例子2

对 sin(x) 做Newton插值逼近

%%
clc;clear all;
x0 = linspace(0,15,15);
y0 = sin(x0);
x= linspace(0,15,100);
y = sin(x);
[s,M] = Newton_interp(x0,y0,x);
plot(x,y,'r');
hold on
    plot(x,s,'b');
    plot(x0,y0,'o');
    hold off
    legend('sin(x)','N(x)');

运行后得到

最后小结

这次写的Newton插值代码只强调了差商计算表对计算的一些简便之处，对于一组相同的样本点$x0,y0$，不同的近似点x,有一张共用的差商计算表，从而节省一定的计算量。

如果是在原有样本点的基础上添加新的样本点，需要注意将其的横坐标按升序重新排列，再结合旧的差商表，节省计算量的同时计算出新的差商表，这里不再进行探讨，如果有兴趣的话，可以自己实现一下。