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
测试依旧用原来的函数生成的数据,即下表:
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.]]
显然这和之前的多项式是一致的。