Python数值计算（21）——非扭结点三次样条曲线

news2025/1/24 5:27:02

前面介绍到紧固和自然三次样条曲线，这次介绍一下非扭结点三次样条曲线。所谓的非扭结点，是指由于最开始的两个子区间使用插值多项式相同，最后两个子区间所使用的插值多项式也相同，这就会导致 $x_1,x_{n-1}$ 在这段多项式上起不到扭结点的效果，故而得名not-a-knot。

本篇除了介绍非扭结点的三次样条曲线外，还将三种不同边界条件的三次样条进行了对比分析。

1. 数学知识

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

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

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 kcsIntp:
    __coef=None
    __bpx=None
    __bpy=None
    def __init__(self,x:np.ndarray,y:np.ndarray):
        '''
        非扭结点三次样条曲线
        '''
        n,=x.shape
        h=np.diff(x)
        a=y.copy()
        dy=np.diff(y)
        A=np.zeros((n,n))         
        A[0,0:3]=[-h[1],h[0]+h[1],-h[0]]
        A[-1,-3:]=[-h[n-2],h[n-3]+h[n-2],-h[n-3]]
        B=np.zeros(n)
        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. 算法测试

采用非扭结点三次样条插值，在 $[0,4]$ 上，对函数 $f(x)=e^x$ 进行拟合，假设我们知道 $x=0,1,2,3,4$ 点处的函数值，以及在 $x=0$ 和 $x=4$ 时的导数值，绘制原函数曲线，以及拟合后的曲线，代码如下：

x=np.array([0,1,2,3,4])
y=np.exp(x)
 
X=np.linspace(0,4,100)
Y=np.exp(X)
plt.plot(X,Y,'r')
 
S=kcsIntp(x,y)
Y1=S(X)
plt.plot(X,Y1,'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，该参数决定了边界条件，使用'not-a-knot'值时，就是非扭结点三次样条曲线，这个参数也是默认的，所实现的拟合效果与第3节中相同，可以通过系数对比确认这一点：

print(sp.c)
'''
[[ 0.48231853  0.48231853  2.66162144  2.66162144]
 [ 0.02929062  1.47624622  2.92320182 10.90806614]
 [ 1.20667268  2.71220952  7.11165756 20.94292553]
 [ 1.          2.71828183  7.3890561  20.08553692]]
'''
print(S.c)
'''
[[ 1.          1.20667268  0.02929062  0.48231853]
 [ 2.71828183  2.71220952  1.47624622  0.48231853]
 [ 7.3890561   7.11165756  2.92320182  2.66162144]
 [20.08553692 20.94292553 10.90806614  2.66162144]]
'''

同样的，系数的组织方式不同，具体细节见前两个章节。

5. 不同类型三次样条曲线的对比

前面根据三种不同边界条件，以下将三种边界条件所构造的分段多项式，对比其1阶，二阶和三阶导数。实现代码如下：

import numpy as np
import matplotlib.pyplot as plt
from numpy.polynomial import Polynomial
from numpy import linalg
from scipy.interpolate import CubicSpline

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

spc=CubicSpline(x,y,bc_type=((1, 1), (1, np.exp(4))))
spn=CubicSpline(x,y,bc_type='natural')
spk=CubicSpline(x,y)

fig,ax=plt.subplots(2,2)
X=np.linspace(0,4,100)
Yc0=spc(X)
Yn0=spn(X)
Yk0=spk(X)
# 拟合函数
ax[0,0].plot(X,Yc0,'r')
ax[0,0].plot(X,Yn0,'g')
ax[0,0].plot(X,Yk0,'b')
ax[0,0].set_title('order-0')
ax[0,0].grid()

# 拟合函数一阶导
Yc1=spc.derivative(1)(X)
Yn1=spn.derivative(1)(X)
Yk1=spk.derivative(1)(X)
ax[0,1].plot(X,Yc1,'r')
ax[0,1].plot(X,Yn1,'g')
ax[0,1].plot(X,Yk1,'b')
ax[0,1].set_title('order-1')
ax[0,1].grid()

# 拟合函数二阶导
Yc2=spc.derivative(2)(X)
Yn2=spn.derivative(2)(X)
Yk2=spk.derivative(2)(X)
ax[1,0].plot(X,Yc2,'r')
ax[1,0].plot(X,Yn2,'g')
ax[1,0].plot(X,Yk2,'b')
ax[1,0].set_title('order-2')
ax[1,0].grid()

# 拟合函数三阶导
Yc3=spc.derivative(3)(X)
Yn3=spn.derivative(3)(X)
Yk3=spk.derivative(3)(X)
ax[1,1].plot(X,Yc3,'r')
ax[1,1].plot(X,Yn3,'g')
ax[1,1].plot(X,Yk3,'b')
ax[1,1].set_title('order-3')
ax[1,1].grid()
plt.show()

对比图如下：