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

效果如下：

可以看到，一阶导（红色）是光滑的，二阶导是连续的（红色），三阶导（蓝色）是分段保持恒定的。