函数插值法之牛顿插值法 python代码实现

老规矩，数学原理什么的就不写了。

直接贴代码和实例演示，以下代码基于python和numpy。

这里比较详细的讲了原理，在这我只关注如何代码实现。

往期博客：

线性方程组的迭代法 python代码实现

数值积分 python代码实现

数值微分 python代码实现

python计算信息熵、信息增益和信息增益率

首先，什么是差商表？就用一张图解释。

差商表是求近似值的前一步，那这里实现求近似值吗？

答案是：现在不。 用代码实现半自动求解。

因为求近似值还要选择节点、代入迭代公式，虽然代码可以实现，但博主觉得“选择节点、代入迭代公式”这个要用代码实现，好像特别麻烦的样子，因为涉及决策的问题，所以还不如自己对照着代入更方便一点 。

如果哪天突然觉得后者更麻烦，不排除用代码实现近似值的求解。
不用哪天，当天当晚就觉得更麻烦。

2020.3.28 23:09

真香！

开始写数学作业后，发现“带入迭代公式”也好麻烦，于是还是用了电脑半自动求解，具体看“求解近似值”。

构建差商表

首先

import numpy as np

定义函数

以下便是我定义的函数：

def csb(x,f,j):
    f0=np.zeros((j+1,x.shape[0]))
    if type(f) is np.ndarray:
        f0[0]=f.copy()
    else:
        for i in range(x.shape[0]):
            f0[0,i]=f(x[i])
    for i in range(1,j+1):
        for k in range(i,j+1):
            f0[i,k]=(f0[i-1,k]-f0[i-1,k-1])/(x[k]-x[k-i])
    f1=np.vstack([x,f0])
    print('所求%d阶差商表如下所示' % j,'\n',f1.T)
    return f1.T

参数说明

“x”是上图中第一列“xk”的值组成的数组。

“f”是上图中的第二列，输入的可以是确定值，也可以是函数，后面会具体演示两者有什么区别。

“j”是最大阶数，在上图中，最大阶数是3阶，那么“j”就是3。

实例运行

拿个实例运行一下，在这个实例中“f”是确定值。

x=np.array([-1,0,1,3])
f=np.array([4,-1,2,6])
csb(x,f,3)

得出以下结果：

所求3阶差商表如下所示
[[-1. 4. 0. 0. 0. ]
[ 0. -1. -5. 0. 0. ]
[ 1. 2. 3. 4. 0. ]
[ 3. 6. 2. -0.33333333 -1.08333333]]

这个数组对照着下图看就知道每一列数字代表着什么意思。

在上个实例中，“f”是确定值。下面再运行一个实例，其中的“f”是函数。

x=np.array([-1,0,1,3])
f=lambda x:x**2
csb(x,f,3)

得出以下结果：

所求3阶差商表如下所示
[[-1. 1. 0. 0. 0.]
[ 0. 0. -1. 0. 0.]
[ 1. 1. 1. 1. 0.]
[ 3. 9. 4. 1. 0.]]

求解近似值

待补充。
补充是不可能补充的，最近是不可能补充的，等忙完这一阵吧。

2020.3.28 23:13

没想到打脸来的这么快。。。

在这里我用了电脑半自动求解，怎么个半自动法呢。

就是将“选择节点、代入迭代公式”中的前半部分自己做了，后半部分交给电脑。

以后有空再写个全自动的。

实例运行

举个例子吧，题目如下。

具体操作如下。

x=np.array([77,78,79,81,82])
y=np.array([2.90256,2.97857,3.06173,3.25530,3.36987])
z=csb(x,y,4)

得到如下结果：

所求4阶差商表如下所示
[[7.70000000e+01 2.90256000e+00 0.00000000e+00 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[7.80000000e+01 2.97857000e+00 7.60100000e-02 0.00000000e+00
0.00000000e+00 0.00000000e+00]
[7.90000000e+01 3.06173000e+00 8.31600000e-02 3.57500000e-03
0.00000000e+00 0.00000000e+00]
[8.10000000e+01 3.25530000e+00 9.67850000e-02 4.54166667e-03
2.41666667e-04 0.00000000e+00]
[8.20000000e+01 3.36987000e+00 1.14570000e-01 5.92833333e-03
3.46666667e-04 2.10000000e-05]]

密密麻麻，看都不好看，那么第一列和第二列不看，就看我需要的那段。

print(z[:,2:6])

[[0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[7.60100000e-02 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[8.31600000e-02 3.57500000e-03 0.00000000e+00 0.00000000e+00]
[9.67850000e-02 4.54166667e-03 2.41666667e-04 0.00000000e+00]
[1.14570000e-01 5.92833333e-03 3.46666667e-04 2.10000000e-05]]

这样好一些了，然后第一问就解决了。

关键是第二问，“取x0=77…x4=82，x=81.12”，然后半自动代码如下。

n=z[0,1]
x1=81.12
for i in range(4):
    s=1
    k=0
    while k<i+1:
        s=s*(x1-x[k])
        k+=1
    n=n+z[i+1,i+2]*s
print(n)

3.2683300909465607

“n”是我们所要求的值，四次的话从头到尾都要。但如果是三次呢？那么我们要刨去第一个值“x=77”，这时的“x0”是“78”，其它的依次递进，然后修改后的半自动代码如下。

n=z[1,1]
x1=81.12
for i in range(3):
    s=1
    k=0
    while k<i+1:
        s=s*(x1-x[k+1])
        k+=1
    n=n+z[i+2,i+2]*s
print(n)

3.268344759040001

我发现《工程与科学计算》这门课，讲真的不难，就是特别特别特别麻烦。