1

我正在尝试对由多个 Bessel 函数（第一类和第二类）组成的实值被积函数进行数值评估。被积函数是振荡衰减的，需要在 0 到 +∞ 之间求值。到目前为止，我使用 scipy.integrate 子包（quad 和 fixed_quad）的尝试都没有成功。当实际上它应该是平滑的时，评估值会跳跃。对于某些参数值集，我还会收到警告：“IntegrationWarning：积分可能发散或缓慢收敛。” （已知是收敛的）或“IntegrationWarning：已达到最大细分数（50）。”

方程式来自：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()
4

1 回答 1

6

我要告诉你一个坏消息：振荡被积函数必须通过特殊方法处理，例如 Filon 方法。阅读 scipy 文档，我认为没有办法做到这一点。但是，似乎MPMath精确地实现了您正在寻找的东西，甚至给出了使用 Bessel 函数的示例。

于 2018-01-20T18:34:22.630 回答