Python数值计算（15）——Newton插值

1. 数学知识

虽然在给定了N个点以后，通过这个点的最小幂多项式是确定的，但是表达方式可不止一种，理论上来说是有无数种的，只要换一组基函数，就可以有一种表示方式。当然最直观常见的，还是前面提到的系数方式，根方式，还有插值的Lagrange形式等。这里介绍另外一种表达方式：

显然这个式子最高次幂不会超过n，然后在时，有：

使用这种表达方式的关键自然是确定前面的系数，而Newton差商公式，就是解决这一问题的。

求解差商的秘诀还是构造差商表，和前面Neville方法类似，但是最后的构造形式会有些差别。

首先，差商的计算方式如下：

零阶差商用函数值表示，最后所求的系数在差商矩阵的主对角线上。

2. 代码实现

Python实现代码如下：

import numpy as np
import matplotlib.pyplot as plt
from numpy.polynomial import Polynomial as P
np.polynomial.set_default_printstyle("unicode")

def newtonIntp(x:np.ndarray,y:np.ndarray):
    '''
    返回基于x,y点对的牛顿插值多项式，以及差商矩阵
    '''
    n=len(x)
    F=np.array([0.0]*n*n).reshape(n,n)
    F[:,0]=y.copy()
    for c in range(1,n):
        for r in range(c,n):
            F[r,c]=(F[r,c-1]-F[r-1,c-1])/(x[r]-x[r-c])
    diag = np.diagonal(F)
    p=P(F[0,0])
    for i in range(1,n):
        p+=F[i,i]*P.fromroots(x[0:i])
    return p,F

3. 测试

测试依旧用原来的函数生成的数据，即下表：

X	Y
-2	-19
-1	0
0	1
1	2
2	21

 测试拟合代码如下：

x=[-2,-1,0,1,2]
y=[-19,0,1,2,21]
Q,M=newtonIntp(x,y)
print(Q)
print(M)

输出为：

1.0 - 2.0·x + 0.0·x² + 3.0·x³
[[-19.   0.   0.   0.   0.]
 [  0.  19.   0.   0.   0.]
 [  1.   1.  -9.   0.   0.]
 [  2.   1.   0.   3.   0.]
 [ 21.  19.   9.   3.   0.]]

显然这和之前的多项式是一致的，实际上可以证明，在过这些点的最小多项式的最简形式（例如升幂排列），有且只有一个，当然，表示方式可以千变万化。