前面介绍到紧固三次样条曲线，这次介绍一下自然三次样条曲线。
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 ncsIntp:
__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]=1
A[-1,-1]=1
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.__bpx3. 算法测试
采用自然三次样条插值，在上，对函数
进行拟合，假设我们知道
点处的函数值，以及在
和
时的导数值，绘制原函数曲线，以及拟合后的曲线，代码如下：
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=ncsIntp(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，该参数决定了边界条件，使用'natural'值时，就是自然三次样条曲线，所实现的拟合效果与第3节中相同，可以通过系数对比确认这一点：
print(sp.c)
'''
[[ 6.07155227e-01 -8.32836920e-02 4.79919365e+00 -5.32306519e+00]
[-2.22044605e-16 1.82146568e+00 1.57161460e+00 1.59691956e+01]
[ 1.11112660e+00 2.93259228e+00 6.32567257e+00 2.38664827e+01]
[ 1.00000000e+00 2.71828183e+00 7.38905610e+00 2.00855369e+01]]
'''
print(S.c)
'''
[[ 1. 1.1111266 0. 0.60715523]
[ 2.71828183 2.93259228 1.82146568 -0.08328369]
[ 7.3890561 6.32567257 1.5716146 4.79919365]
[20.08553692 23.86648273 15.96919556 -5.32306519]]
'''同样的，系数的组织方式不同，具体细节见前1章节。
