前面介绍到紧固和自然三次样条曲线,这次介绍一下非扭结点三次样条曲线。所谓的非扭结点,是指由于最开始的两个子区间使用插值多项式相同,最后两个子区间所使用的插值多项式也相同,这就会导致在这段多项式上起不到扭结点的效果,故而得名not-a-knot。
本篇除了介绍非扭结点的三次样条曲线外,还将三种不同边界条件的三次样条进行了对比分析。
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 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. 算法测试
采用非扭结点三次样条插值,在上,对函数
进行拟合,假设我们知道
点处的函数值,以及在
和
时的导数值,绘制原函数曲线,以及拟合后的曲线,代码如下:
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()
对比图如下:
上图中,红色、绿色和蓝色分别表示的是紧固、自然和非扭结点三次样条曲线,order为0~3分别表示0阶,1阶,2阶和3阶导数。可以看到:
1. 拟合曲线在一阶导数上仍旧是光滑的,在二阶导数上是连续的,但是在三阶导数上就存在跳跃了,不过,在每个子区间内,其导函数值保持恒定,这些都是三次样条曲线应该具有的基本特征。
2. 对于非扭结点三次样条曲线,可以看到在右下图中,前两个子区间的三阶导数值是相同的,而另外两种三次样条曲线则不具备这一特性。
3. 由于紧固三次样条曲线使用了更多的原函数信息,因此,其1阶,2阶的趋势基本上和原函数保持一致,而自然三次样条曲线在二阶导数中的表现是最差的,在三阶导数中也产生了较大的振动。