1. 概述
不管是前面介绍到拉格朗日插值还是牛顿插值，拟合的函数比线性插值更加“优秀”，即它们都是连续可导的，但是，有时拟合还有这样的要求，就是除了在给定点处的函数值要相等外，还要求在这些指定点处的导数值也要是指定的值，甚至是m阶导数相等，也就是所谓的密切多项式，处理这类拟合的插值方法，就是Hermite插值。
2. 数学原理
如果有一个位置函数，和严格递增的序列
在区间
上，在点
与
，
都相同，则次数最小的多项式是唯一的，被称为Hermite多项式，其次数最多为
，并由下面的式子给出：
其中，若用表示第
个
次的Lagrange插值多项式，则有：
3. 算法实现
根据前述的数学原理，实现Hermite插值算法如下：
import numpy as np
import matplotlib.pyplot as plt
from numpy.polynomial import Polynomial
np.polynomial.set_default_printstyle("unicode")
def hermiteIntp(x,y,dy):
'''
返回Hermite插值多项式
'''
n=len(x)
H=Polynomial([0])
for k in range(n):
roots=list(x)
s=roots[k]
del roots[k]
L=Polynomial.fromroots(roots)
L/=L(s)
dL=L.deriv()
p1=Polynomial([-s,1])
L2=L**2
H+=y[k]*(1-2*p1*dL(s))*L2
H+=dy[k]*p1*L2
return H
x=np.array([-1,0,1])
y=np.array([9,3,1])
dy=np.array([-6,-2,-6])
H=hermiteIntp(x,y,dy)
X=np.linspace(-1,1,50)
Y=H(X)
plt.plot(X,Y,'r')
# 以上拟合的数据由如下多项式产生
a=Polynomial([3,-2,4,-3,-2,1])
plt.plot(X,a(X),'g.')
plt.grid()
plt.show()显示图形如下：
可见两者是完全贴合的。
不过，由于原函数本来就是一个多项式，这没什么好奇怪的，改用非多项式函数进行拟合测试，例如函数
，手工计算其导函数为
，编写代码如下：
#原函数
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) # 一阶导数值
H=hermiteIntp(x,y,dy)
X=np.linspace(0,1,50)
Y=H(X)
plt.plot(X,Y,'r')
# 绘制原函数
plt.plot(X,f(X),'b-.')
plt.grid()
plt.show()绘制图形如下：
可以看到两者的吻合度还是很高的。
4. 使用库
SciPy里的Hermite插值实现函数是scipy.interpolate模块下的KroghInterpolator类，其构造函数X和Y必须同等长度，可以给出若干个坐标点作为形参，也可以给出每个x对应的函数和导数值，KroghInterpolator要求x是递增序列，如果X中多次某一个x轴坐标xi出现多次，例如k次，那么与之对应的y分别是该点处的函数值，一阶导函数值，……（k-1）阶导函数值，对于Hermite插值而言，就是每个X中的元素写两次，前面的例子可以用如下代码实现：
import numpy as np
from numpy.polynomial import Polynomial
import matplotlib.pyplot as plt
from scipy.interpolate import KroghInterpolator
x=np.array([-1,-1,0,0,1,1])
y=np.array([9,-6,3,-2,1,-6])
H=KroghInterpolator(x,y)
X=np.linspace(-1,1,50)
Y=H(X)
plt.plot(X,Y,'r')
# 以上拟合的数据由如下多项式产生
a=Polynomial([3,-2,4,-3,-2,1])
plt.plot(X,a(X),'g.')
plt.grid()
plt.show()运行结果显示的图形和原来相同。
