Python数值计算（17）——Hermite插值

这次介绍一下使用差商表构造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()

结果如下图所示：

这个图形和前面一节所介绍的方法完全相同，可见该算法是有效的。