Python数值计算（19）——紧固三次样条曲线

前面介绍到了三种样条曲线的类型，接下来将分三个章节分别介绍紧固三次样条曲线，自然三次样条曲线和非扭结点三次样条曲线。

1. 数值知识

这个在前面已经做过介绍，这里再次重复说明一遍，它对我们算法实现具有很大的帮助：

同样的，就是各分段起点的函数值，通过上述方程组求解出每一个以后，可以把计算出来。

2. 算法实现

和线性插值一样，我们用一个类来实现插值算法，以及其他的方法、属性，这样就可以在后续实现方便的多值估算和绘图，最终实现代码如下：

import numpy as np
import matplotlib.pyplot as plt
from numpy.polynomial import Polynomial
from numpy import linalg
from scipy.interpolate import CubicSpline
np.polynomial.set_default_printstyle("unicode")

class ccsIntp:
    __coef=None
    __bpx=None
    __bpy=None
    def __init__(self,x:np.ndarray,y:np.ndarray,fp0,fpn):
        '''
        紧固三次样条曲线,fp0和fpn是端点处的导数值
        '''
        n,=x.shape
        h=np.diff(x)
        a=y.copy()
        dy=np.diff(y)
        A=np.zeros((n,n)) 
        A[0,0:2]=[2*h[0],h[0]]
        A[-1,-2:]=[h[-2],2*h[-2]]
        B=np.zeros(n)
        B[0]=3*dy[0]/h[0]-3*fp0
        B[n-1]=3*fpn-3*dy[n-2]/h[n-2]
        for i in range(1,n-1):
            A[i,i-1:i+2]=[h[i-1],2*(h[i-1]+h[i]),h[i]]
            B[i]=3*dy[i]/h[i]-3*dy[i-1]/h[i-1]
        c=linalg.solve(A,B)
        d=np.zeros(n)
        b=np.zeros(n)
        for j in range(n-1):
            d[j]=(c[j+1]-c[j])/3/h[j]
            b[j]=dy[j]/h[j]-h[j]/3*(2*c[j]+c[j+1])
        self.__coef=np.array([a,b,c,d])[:,:-1].T
        self.__bpx=x.copy()
        self.__bpy=y.copy()

    def __call__(self,w):
        n,=w.shape
        ret= np.zeros(n)
        for i in range(n):
            if self.__bpx[0]>=w[i]:
                ret[i]=self.__bpy[0]
                continue
            if self.__bpx[-1]<=w[i]:
                ret[i]=self.__bpy[-1]
                continue
            j=0
            while self.__bpx[j]<w[i]:
                j+=1
            cp=self.__coef[j-1,:]
            p=Polynomial([0])
            for t in range(len(cp)):
                p+=cp[t]*Polynomial([-self.__bpx[j-1],1])**t
            ret[i]=p(w[i])
        return ret
    
    @property
    def c(self):
        '''
        如果提供的是n+1个点对，则系数是shape为(n,4)的ndarray
        每一行就是一个分段区间的参数，依次记作a,b,c,d
        则该区间的样条曲线就是y=a+b*(x-xj)+c*(x-xj)**2+d*(x-xj)**3
        其中0<=j<=n-1    
        '''
        return self.__coef
    @property
    def x(self):
        return self.__bpx

3. 算法测试

采用紧固三次样条插值，在上，对函数进行拟合，假设我们知道点处的函数值，以及在和时的导数值，绘制原函数曲线，以及拟合后的曲线，代码如下：

x=np.array([0,1,2,3,4])
y=np.exp(x)

S=ccsIntp(x,y,1,np.exp(4))# 后两个值是e^x在0和4处的导数值

# 绘制图形
X=np.linspace(0,4,100)
Y1=np.exp(X)
plt.plot(X,Y1,'r')
Y2=S(X)
plt.plot(X,Y2,'b-.')
plt.grid()
plt.show()

运行效果如下：

可以看到，这两者的图形几乎完全贴合在一起。

4. 现有工具包

在scipy的工具包中，scipy.interpolate.CubicSpline类可以完成三次样条曲线的插值功能，构造函数原型如下：

class CubicSpline(x, y, axis=0, bc_type='not-a-knot', extrapolate=None)

其中x，y是拟合点，主要的区别是bc_type，该参数决定了边界条件，使用'clamped'值时，就是紧固三次样条曲线，但是，需要注意的是，在Scipy中，紧固边界条件是在端点处的一阶导是为0，而不是向上面那样的使用给定值，如果需要使用给定值，则应该使用如下的方式：

bc_type=((1, fp0), (1, fpn))

其中fp0和fpn就是端点处的导数值，测试代码如下：

x=np.array([0,1,2,3,4])
y=np.exp(x)
sp=CubicSpline(x,y,bc_type=((1, 1), (1, np.exp(4))))
S=ccsIntp(x,y,1,np.exp(4))
X=np.linspace(0,4,100)
Y1=np.exp(X)
plt.plot(X,Y1,'r')
Y3=sp(X)
plt.plot(X,Y3,'b-.')
plt.grid()
plt.show()

其绘制图形效果和我们前面自定义类实现的效果相同。

我们还可以通过其属性了解分段拟合多项式的系数，对比我们自定义函数和CubicSpline实例的系数：

x=np.array([0,1,2,3,4])
y=np.exp(x)
sp=CubicSpline(x,y,bc_type=((1, 1), (1, np.exp(4))))
S=ccsIntp(x,y,1,np.exp(4))
print(sp.c)
'''
[[ 0.26217911  0.72939145  1.89346921  5.48715745]
 [ 0.45610272  1.24264005  3.4308144   9.11122202]
 [ 1.          2.69874277  7.37219722 19.91423364]
 [ 1.          2.71828183  7.3890561  20.08553692]]
'''
print(S.c)
'''
[[ 1.          1.          0.45610272  0.26217911]
 [ 2.71828183  2.69874277  1.24264005  0.72939145]
 [ 7.3890561   7.37219722  3.4308144   1.89346921]
 [20.08553692 19.91423364  9.11122202  5.48715745]]
'''

可以看到，我们自定义的类，系数是按行升幂排列的，即第行就是分段的系数，而CubicSpline类的系数，则是按列降幂排列的即第列就是分段的系数，因此使用时需要注意区分，无论如何，通过两种方式，都可以很快的写出分段多项式方程：