Python数值计算（24）——PCHIP
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 pltdef 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**2H+=y[k]*(1-2*p1*dL(s))*L2H+=dy[k]*p1*L2return 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=0elif (np.sign(d1)!=np.sign(d2)) and (np.abs(d)>np.abs(2*d1)):d=3*d1return ddef __init__(self,x,y):self.__bpx=x.copy()self.__bpy=y.copy()n=len(x)h=np.diff(x)d=np.diff(y)/hself.__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])continueif w[i]>=self.__bpx[-1]:ret[i]=self.__polys[-1](w[i])continuej=0while w[i] > self.__bpx[j]:j+=1ret[i]=self.__polys[j-1](w[i])return ret@propertydef x(self):return self.__bpx@propertydef y(self):return self.__bpy@propertydef dy(self):return self.__dy@propertydef c(self):coef=np.zeros((len(self.__polys),4))for i in range(len(self.__polys)):coef[i,:]=self.__polys[i].coefreturn coef3. 算法测试
可以测试前面生成的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()效果如下：
可以看到，一阶导（红色）是光滑的，二阶导是连续的（红色），三阶导（蓝色）是分段保持恒定的。