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

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()

运行结果显示的图形和原来相同。