1. 算法描述
PCHIP是“分段三次埃尔米特插值多项式(piecewise cubic Hermite interpolating polynomial)”的缩写,但是单纯根据这个名字,其实并不能确定是哪一种分段三次埃尔米特插值多项式,因为样条插值也是分段三次埃尔米特插值多项式的一种,只是对斜率的约束条件不同而已,这里讨论的是一种保形(shape preserving)的特定插值函数,和Akima插值算法类似,其关键的思想还是如何确定各个采样点的斜率,使函数的值不会再局部上过多的偏离给定的数据。如果该点和前面一个点的差分值符号相反,或者有一个为零,则在
处函数为离散的极大值或者极小值,因此可以令
,如果
符号相同,则可以令:
但是对于首位两个端点处的斜率,需要用到一点不同的算法,详细的实现见算法描述部分。
2. 算法实现
首先,为了构造Hermite插值多项式,我们要用到前面实现的Hermite插值函数:
import numpy as np
from numpy.polynomial import Polynomial
from scipy.interpolate import PchipInterpolator,CubicHermiteSpline
import matplotlib.pyplot as plt
def hermiteIntp(x,y,dy):
'''
通过Lagrange插值法返回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
然后,我们可以定义一个pchipIntp的类,通过它保持分段插值多项式的信息,并为后面计算其他点的函数值提供便利。注意看其中包含了一个私有的成员方法, __enddf,专门用来计算首位端点处的斜率。
class pchipIntp:
__bpx:np.ndarray # 插值点
__bpy:np.ndarray # 插值点处函数值
__dy:np.ndarray # 插值点处导数值
__polys=[] # 分段多项式
def __enddf(self,h1,h2,d1,d2):
'''
计算两个端点处的斜率值
'''
d=((2*h1+h2)*d1-h1*d2)/(h1+h2)
if np.sign(d) != np.sign(d1):
d=0
elif (np.sign(d1)!=np.sign(d2)) and (np.abs(d)>np.abs(2*d1)):
d=3*d1
return d
def __init__(self,x,y):
self.__bpx=x.copy()
self.__bpy=y.copy()
n=len(x)
h=np.diff(x)
d=np.diff(y)/h
self.__dy=np.zeros(n)
# 首位处的导数值
self.__dy[0]=self.__enddf(h[0],h[1],d[0],d[1])
self.__dy[-1]=self.__enddf(h[-1],h[-2],d[-1],d[-2])
# 中间各点导数值
for k in range(1,n-1):
if d[k-1]*d[k]>0:
w1=2*h[k]+h[k-1]
w2=h[k]+2*h[k-1]
self.__dy[k]=(w1+w2)/(w1/d[k-1]+w2/d[k])
# 构造分段Hermite多项式,使用了前面的算法
for k in range(n-1):
# 每个子区间使用一个首位端点的函数值和导数值,
# 因此该区间段的最高次数为3次
p=hermiteIntp(self.__bpx[k:k+2],self.__bpy[k:k+2],self.__dy[k:k+2])
self.__polys.append(p)
def __call__(self,w:np.ndarray):
n=len(w)
ret=np.zeros(n)
for i in range(n):
if w[i]<=self.__bpx[0]:
ret[i]=self.__polys[0](w[i])
continue
if w[i]>=self.__bpx[-1]:
ret[i]=self.__polys[-1](w[i])
continue
j=0
while w[i] > self.__bpx[j]:
j+=1
ret[i]=self.__polys[j-1](w[i])
return ret
@property
def x(self):
return self.__bpx
@property
def y(self):
return self.__bpy
@property
def dy(self):
return self.__dy
@property
def c(self):
coef=np.zeros((len(self.__polys),4))
for i in range(len(self.__polys)):
coef[i,:]=self.__polys[i].coef
return coef
3. 算法测试
可以测试前面生成的pchip插值曲线:
x=np.array([0,1,2,3,4,5,6])
y=np.array([5,4,0,4,6,1,2])
xp=pchipIntp(x,y)
print(xp.c)
t=np.linspace(x[0]-1,x[-1]+1,50)
plt.plot(t,xp(t),'b-.')
输出的系数如下:
[[ 5.00000000e+00 0.00000000e+00 -1.40000000e+00 4.00000000e-01]
[-9.60000000e+00 3.52000000e+01 -2.80000000e+01 6.40000000e+00]
[ 8.00000000e+01 -1.01333333e+02 4.13333333e+01 -5.33333333e+00]
[ 3.80000000e+01 -3.73333333e+01 1.26666667e+01 -1.33333333e+00]
[-8.74000000e+02 6.00000000e+02 -1.35000000e+02 1.00000000e+01]
[-1.24000000e+02 7.50000000e+01 -1.50000000e+01 1.00000000e+00]]
每一行表示一个三次插值多项式,注意这里没有考虑起点的问题,因此,任何一段之间的函数都是:
绘制效果如下:
4. 现有工具包
Scipy中其实有两个类都支持这种pchip插值,分别是PchipInterpolator和CubicHermiteSpline,但是具体使用方式上存在差别,前者只需要提供即可完成插值工作,而后者需要第三个参数,亦即各点处的斜率,结合前面的自定义类,测试如下:
x=np.array([0,1,2,3,4,5,6])
y=np.array([5,4,0,4,6,1,2])
xp=pchipIntp(x,y)
t=np.linspace(x[0]-1,x[-1]+1,50)
plt.plot(t,xp(t),'b-.')
p=PchipInterpolator(x,y)
plt.plot(t,p(t),'r--')
chs=CubicHermiteSpline(x,y,xp.dy)
print(p.c-chs.c)
plt.plot(t,chs(t),'c.')
plt.grid()
plt.show()
绘制效果如下:
可以看到三者是一致的。并且输出的系数差:
[[0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0.]]
可见两者的系数是完全相同的。 但是它们的系数都是采用分段方式,而且是按幂升幂排列的,即有:
另外,还可以查看一下插值后的1,2,3阶导函数情况:
fig=plt.figure()
plt.plot(t,p.derivative(1)(t),'r')
plt.plot(t,p.derivative(2)(t),'g')
plt.plot(t,p.derivative(3)(t),'b.')
plt.grid()
plt.show()
效果如下:
可以看到,一阶导(红色)是光滑的,二阶导是连续的(红色),三阶导(蓝色)是分段保持恒定的。