我正在进行数值积分,其中要积分的函数使用三次样条表示。三次样条曲线在函数MEcompute中作为splc启动

现在实际执行插值的integrand需要三次样条数组,因此我需要将splc传递给此新函数。在这里我被困住了。

# function defining the integrand which uses the spline coef array to give interpolated values
def integrand(xpoint):
    spline_array=splc
    result=interpolate.splev(xpoint,spline_array,der=0)
    return result

#----------------------------------------
# function to the matrix element for certain rovibrational state
def MEcompute(psi1,psi2,psi_r, parameter, parameter_r ):

        # step 1: gen cubic spline coefs.
        splc=interpolate.splrep(parameter_r,parameter,s=0)

        # generate interpolated parameter for same xaxis as psi
        parameter_interp=interpolate.splev(psi_r,splc,der=0)

        # compute the pointwise products
        p1=np.multiply(psi1,psi2)
        p2=np.multiply(p1,psi_r)
        p3=np.multiply(p2,psi_r)
        product=np.multiply(p3,parameter_interp)

        # step 1: gen cubic spline coefs
        splc=interpolate.splrep(psi_r,product,s=0)

        # compute the integral using adaptive Quadrature
        #result=integrate.quadrature(integrand,0.2,4.48,tol=1.0e-9,maxiter=500)
        result=integrate.quadrature(integrand,0.2,4.48,tol=1.0e-9,maxiter=500)
        print("<psi1|parameter|psi2> = ",result)

#----------------------------------------
# computing the value

MEcompute(v1,v2,rwave,parameter1,distance)
#----------------------------------------


我收到错误,

NameError: name 'splc' is not defined


因为integrand函数看不到在函数splc中启动的MEcompute数组,所以发生了这种情况。

现在我有一个克服这个想法的想法:


splc中将数组MEcompute导出为txt文件,然后在integrand函数中加载该txt文件。这肯定会增加计算时间。


有人可以建议一种更好的方法来做到这一点。

最佳答案

使用args= keyword argument将额外的参数传递给函数以进行集成:

result = integrate.quadrature(integrand, 0.2, 4.48,
                              tol=1.0e-9, maxiter=500,
                              args=(splc,))


并修改您的被积以接受参数:

def integrand(xpoint, splc):
    spline_array=splc
    result=interpolate.splev(xpoint,spline_array,der=0)
    return result

关于python - 将函数内生成的数组传递给python中另一个称为函数(集成)的对象,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/48062522/

10-14 19:03
查看更多