Python数值计算（20）——自然三次样条曲线

C囧囧

已于 2024-08-04 14:07:54 修改

阅读量453

点赞数 8

分类专栏： Python 文章标签： python 开发语言

于 2024-08-04 13:40:53 首次发布

本文链接：https://blog.csdn.net/cdinten/article/details/140893421

版权

Python 专栏收录该内容

28 篇文章 0 订阅

订阅专栏

前面介绍到紧固三次样条曲线，这次介绍一下自然三次样条曲线。

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 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.__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=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章节。