python 中高度振荡的一维被积函数（包含贝塞尔函数）的数值积分

Question

我正在尝试对由多个贝塞尔函数（第一类和第二类）组成的实值被积函数进行数值计算。被积函数振荡和衰减，需要在 0 和 +∞ 之间求值。到目前为止，我尝试使用 scipy.integrate 子包（quad 和 fixed_quad）都没有成功。评估值跳来跳去，而实际上它应该是平滑的。对于某些参数值集，我还收到警告："IntegrationWarning: The integral is probably divergent, or slowly convergent."（已知是收敛的）或"IntegrationWarning: The maximum number of subdivisions (50) has been achieved."

等式来自：http://dx.doi.org/10.1029/WR003i001p00241

也可以在这里找到：http://www.aqtesolv.com/papadopu.htm

感谢您在 python...

中对繁琐函数的数值积分提供的任何帮助

代码示例

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import special as sps
import scipy.integrate as integrate

# define constants and variables (SI mks units):
r_w = 0.15
r_c = 0.16
b = 10
S_s = 1E-6
Q = 0.001
S = S_s*b
K=1E-8
T=K*b
alpha = (r_w**2)*S/r_c**2
def r_D(r):
    return r/r_w
def u(r,t):
    return r**2*S/(4*T*t)
def B(beta):
    return beta*sps.jv(0,beta) - 2*alpha*sps.jv(1,beta)
def A(beta):
    return beta*sps.yn(0,beta) - 2*alpha*sps.yn(1,beta)
def intd_f(beta,r,t): 
    TOP = (1-np.exp(-beta**2*r_D(r)**2)/(4*u(r,t)))*( sps.jv(0,beta*r_D(r))*A(beta) - sps.yn(0,beta*r_D(r))*B(beta) )
    BOT = (A(beta)**2 + B(beta)**2)*beta**2
    return TOP/BOT
def s(r,t):
    banana = (2*alpha*Q)/(np.pi**2*K*b)
    apple = integrate.quad(intd_f, 0, np.inf, args=(r,t))[0]
    #apple = integrate.fixed_quad(intd_f, 0, 1E100, args=(r,t))[0] # another option I have tried
    return banana*apple

#%% simple example usage...
r=np.arange(1,10,.1)
t=60*60*24*pd.Series([1/24,1,365,3650])
plt.figure()
for tt in t:
    print('time=',tt)
    snow=[]
    for rr in r:
        print('r=',rr)
        snow.append(s(rr,tt))
    plt.subplot(2,1,1)
    plt.plot(r,snow,label='t='+str(tt/(60*60*24))+'d')
    plt.subplot(2,1,2)
    plt.semilogy(r,np.abs(snow))    
plt.subplot(2,1,1)
plt.legend()

Answer 1

有个坏消息要告诉你：振荡被积函数必须用特殊方法处理，例如 Filon 方法。通读 scipy 文档，我看不出有什么办法可以做到这一点。然而，MPMath 似乎准确地实现了您正在寻找的东西，甚至给出了一个使用贝塞尔函数的示例。

python 中高度振荡的一维被积函数（包含贝塞尔函数）的数值积分

Numerical integration of highly oscillating 1-D integrand (containing Bessel functions) in python

python

scipy

numerical-integration

quad

bessel-functions