python – 从UnivariateSpline对象获取样条方程

概述我正在使用UnivariateSpline为我拥有的某些数据构建分段多项式.然后我想在其他程序中使用这些样条(在C或FORTRAN中),因此我想了解生成样条函数背后的等式. 这是我的代码： import numpy as npimport scipy as spfrom scipy.interpolate import UnivariateSplineimport matplotlib.p 我正在使用UnivariateSpline为我拥有的某些数据构建分段多项式.然后我想在其他程序中使用这些样条(在C或FORTRAN中),因此我想了解生成样条函数背后的等式.

这是我的代码：

import numpy as npimport scipy as spfrom  scipy.interpolate import UnivariateSplineimport matplotlib.pyplot as pltimport bisectdata = np.loadtxt('test_C12H26.dat')TmID = 800.0print "TmID",TmIDnmID = bisect.bisect(data[:,0],TmID)fig = plt.figure()plt.plot(data[:,data[:,7],ls='',marker='o',markevery=20)npts = len(data[:,0])#print "npts",nptsw = np.ones(npts)w[0] = 100w[nmID] = 100w[npts-1] = 100spline1 = UnivariateSpline(data[:nmID,data[:nmID,s=1,w=w[:nmID])coeffs = spline1.get_coeffs()print coeffsprint spline1.get_knots()print spline1.get_resIDual()print coeffs[0] + coeffs[1] * (data[0,0] - data[0,0]) \                + coeffs[2] * (data[0,0])**2 \                + coeffs[3] * (data[0,0])**3,\      data[0,7]print coeffs[0] + coeffs[1] * (data[nmID,0]) \                + coeffs[2] * (data[nmID,0])**2 \                + coeffs[3] * (data[nmID,\      data[nmID,7]print TmID,data[-1,0]spline2 = UnivariateSpline(data[nmID-1:,data[nmID-1:,w=w[nmID-1:])print spline2.get_coeffs()print spline2.get_knots()print spline2.get_resIDual()plt.plot(data[:,spline1(data[:,0]))plt.plot(data[:,spline2(data[:,0]))plt.savefig('test.png')

这是由此产生的情节.我相信每个区间都有有效的样条线,但看起来我的样条方程不正确…我找不到任何关于scipy文档中应该是什么的引用.有人知道吗？谢谢 ！

解决方法 关于如何获取系数并手动生成样条曲线,scipy documentation没有任何说法.但是,有可能从现有的B样条文献中找出如何做到这一点.以下函数bspleval显示了如何构造B样条基函数(代码中的矩阵B),通过将系数与最高阶基函数相乘并求和,可以很容易地生成样条曲线：

def bspleval(x,knots,coeffs,order,deBUG=False):    '''    Evaluate a B-spline at a set of points.    Parameters    ----------    x : List or ndarray        The set of points at which to evaluate the spline.    knots : List or ndarray        The set of knots used to define the spline.    coeffs : List of ndarray        The set of spline coefficIEnts.    order : int        The order of the spline.    Returns    -------    y : ndarray        The value of the spline at each point in x.    '''    k = order    t = knots    m = alen(t)    npts = alen(x)    B = zeros((m-1,k+1,npts))    if deBUG:        print('k=%i,m=%i,npts=%i' % (k,m,npts))        print('t=',t)        print('coeffs=',coeffs)    ## Create the zero-order B-spline basis functions.    for i in range(m-1):        B[i,:] = float64(logical_and(x >= t[i],x < t[i+1]))    if (k == 0):        B[m-2,-1] = 1.0    ## Next iteratively define the higher-order basis functions,working from lower order to higher.    for j in range(1,k+1):        for i in range(m-j-1):            if (t[i+j] - t[i] == 0.0):                first_term = 0.0            else:                first_term = ((x - t[i]) / (t[i+j] - t[i])) * B[i,j-1,:]            if (t[i+j+1] - t[i+1] == 0.0):                second_term = 0.0            else:                second_term = ((t[i+j+1] - x) / (t[i+j+1] - t[i+1])) * B[i+1,:]            B[i,j,:] = first_term + second_term        B[m-j-2,-1] = 1.0    if deBUG:        plt.figure()        for i in range(m-1):            plt.plot(x,B[i,k,:])        plt.Title('B-spline basis functions')    ## Evaluate the spline by multiplying the coefficIEnts with the highest-order basis functions.    y = zeros(npts)    for i in range(m-k-1):        y += coeffs[i] * B[i,:]    if deBUG:        plt.figure()        plt.plot(x,y)        plt.Title('spline curve')        plt.show()    return(y)

为了举例说明如何将其与Scipy现有的单变量样条函数一起使用,以下是一个示例脚本.这将获取输入数据并使用Scipy的功能以及面向对象的样条拟合方法.从两者中的任何一个中获取系数和结点,并使用这些作为我们手动计算的例程bspleval的输入,我们重现它们所做的相同曲线.请注意,手动评估曲线与Scipy评估方法之间的差异非常小,几乎可以肯定是浮点噪声.

x = array([-273.0,-176.4,-79.8,16.9,113.5,210.1,306.8,403.4,500.0])y = array([2.25927498e-53,2.56028619e-03,8.64512988e-01,6.27456769e+00,1.73894734e+01,3.29052124e+01,5.14612316e+01,7.20531200e+01,9.40718450e+01])x_nodes = array([-273.0,-263.5,-234.8,-187.1,-120.3,-34.4,70.6,194.6,337.8,500.0])y_nodes = array([2.25927498e-53,3.83520726e-46,8.46685318e-11,6.10568083e-04,1.82380809e-01,2.66344008e+00,1.18164677e+01,3.01811501e+01,5.78812583e+01,9.40718450e+01])## Now get scipy's spline fit.k = 3tck = splrep(x_nodes,y_nodes,k=k,s=0)knots = tck[0]coeffs = tck[1]print('knot points=',knots)print('coefficIEnts=',coeffs)## Now try scipy's object-orIEnted version. The result is exactly the same as "tck": the knots are the## same and the coeffs are the same,they are just querIEd in a different way.uspline = UnivariateSpline(x_nodes,s=0)uspline_knots = uspline.get_knots()uspline_coeffs = uspline.get_coeffs()## Here are scipy's native spline evaluation methods. Again,"ytck" and "y_uspline" are exactly equal.ytck = splev(x,tck)y_uspline = uspline(x)y_knots = uspline(knots)## Now let's try our manually-calculated evaluation function.y_eval = bspleval(x,deBUG=False)plt.plot(x,ytck,label='tck')plt.plot(x,y_uspline,label='uspline')plt.plot(x,y_eval,label='manual')## Next plot the knots and nodes.plt.plot(x_nodes,'ko',markersize=7,label='input nodes')            ## nodesplt.plot(knots,y_knots,'mo',markersize=5,label='tck knots')    ## knotsplt.xlim((-300.0,530.0))plt.legend(loc='best',prop={'size':14})plt.figure()plt.Title('difference')plt.plot(x,ytck-y_uspline,label='tck-uspl')plt.plot(x,ytck-y_eval,label='tck-manual')plt.legend(loc='best',prop={'size':14})plt.show()

总结

以上是内存溢出为你收集整理的python – 从UnivariateSpline对象获取样条方程全部内容，希望文章能够帮你解决python – 从UnivariateSpline对象获取样条方程所遇到的程序开发问题。

如果觉得内存溢出网站内容还不错，欢迎将内存溢出网站推荐给程序员好友。