这次介绍一下使用差商表构造Hermite多项式的方法。
前面介绍到了两种很经典的插值多项式，即Lagrange和Newton插值多项式，并在前一篇中阐述了如何通过Lagrange插值方式构造Hermite多项式，这次通过牛顿差商法构造Hermite多项式。
1. 数学原理
以三个点为例，以及在这3个点处的导数值
，构造差商表如下：
可以看到，差商表的第一列是x坐标值，并且写了两次，第二列则是对应的函数值，也是写了两次，到第三列为一阶差商时（导数值），对于已知点处导数，直接用现有的导数值，而对于没有现成的导数值，则使用类似Newton差分方式计算，随后二阶，三阶等高阶差商依次计算，最后得到主对角线元素值Q，最后如下等式构造Hermite多项式：
推广：对于共
个点，则Hermite多项式如下：
2. 算法实现
根据上面的原理，实现算法如下：
import numpy as np
import matplotlib.pyplot as plt
from numpy.polynomial import Polynomial
def hermiteIntp_dft(x:np.ndarray,y:np.ndarray):
'''
使用差商表计算Hermite插值多项式
'''
n,=x.shape
# differential table
dft=np.zeros((n,n+1))
dft[:,0]=x
dft[0::2,1]=y[0::2]
dft[1::2,1]=y[0::2]
dft[1::2,2]=y[1::2] # 已知导数值
dft[2::2,2]=(dft[2::2,1]-dft[1:-1:2,1])/(dft[2::2,0]-dft[1:-1:2,0]) # 差商值
for c in range(3,n+1):
for r in range(c-1,n):
dft[r,c]=(dft[r,c-1]-dft[r-1,c-1])/(dft[r,0]-dft[r-c+1,0])
Q=np.diag(dft[:,1::])
H=Polynomial([0])
p=Polynomial([1])
for i in range(n):
H+=Q[i]*p
p*=Polynomial([-x[i],1])
return H算法测试
还是以前面的多项式为例进行测试，已知在
时，
，以及对应的导数值
，构造Hermite的代码如下：
x=np.array([-1,-1,0,0,1,1])
y=np.array([9,-6,3,-2,1,-6])
H=hermiteIntp_dft(x,y)
print(H) # 3.0 - 2.0·x + 4.0·x² - 3.0·x³ - 2.0·x⁴ + 1.0·x⁵函数图形此处从略。
再以函数，手工计算其导函数为
，先通过给定点生成函数值和导函数值，然后用这些值进行拟合。并绘制图形比较，代码如下：
#原函数
def f(x):
return np.exp(-2*x)*np.sin(3*x)
#导函数
def df(x):
return np.exp(-2*x)*(-2*np.sin(3*x)+3*np.cos(3*x))
x=np.array([0,0.5,1])
yh=np.vectorize(f)
y=yh(x) # 函数值
dyh=np.vectorize(df)
dy=dyh(x) # 一阶导数值
hx=np.zeros(2*len(x))
hx[0::2]=x
hx[1::2]=x
hy=np.zeros(2*len(x))
hy[0::2]=y
hy[1::2]=dy
# 调用函数，生成插值多项式
H=hermiteIntp_dft(hx,hy)
# 绘制图形
X=np.linspace(0,1,50)
Y=H(X)
plt.plot(X,Y,'r')
# 绘制原函数
plt.plot(X,f(X),'b-.')
plt.grid()
plt.show()
结果如下图所示：
这个图形和前面一节所介绍的方法完全相同，可见该算法是有效的。
