具有无限极限的四定积分

Question

我正在尝试解决以下问题：

​

其中，g是一个常量，\mu = 0，但是我得到了一个警告，结果并不完全是预期的结果。

代码如下：

%matplotlib inline
from scipy.integrate import quad
from numpy import arange
from numpy import pi
import matplotlib.pyplot as plt
import numpy as np

###Parameters

ms = 100. ## m
ggX = 2.  ## g

def Integrand(E, T, m_dm):
    return ((ggX/(2*(np.pi**2)))*E*(E**2-m_dm**2)**(1/2))/(np.exp(E/T)-1)

def Integrate_E(T,m_dm):
    '''
        Integrate over E given T and ms
    '''
    return quad(Integrand, m_dm, np.inf, args=(T,m_dm))[0]

TT = np.logspace(-10,15,100)
nn =[Integrate_E(T,ms) for T in (TT)]

plt.loglog(TT, nn)
plt.grid(True)

这是我得到的结果：

/home/vcti/anaconda3/lib/python3.7/site-packages/scipy/integrate/quadpack.py:385: IntegrationWarning: The algorithm does not converge.  Roundoff error is detected
  in the extrapolation table.  It is assumed that the requested tolerance
  cannot be achieved, and that the returned result (if full_output = 1) is 
  the best which can be obtained.
  warnings.warn(msg, IntegrationWarning)
/home/vcti/anaconda3/lib/python3.7/site-packages/scipy/integrate/quadpack.py:385: IntegrationWarning: The integral is probably divergent, or slowly convergent.
  warnings.warn(msg, IntegrationWarning)
/home/vcti/anaconda3/lib/python3.7/site-packages/scipy/integrate/quadpack.py:385: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
  warnings.warn(msg, IntegrationWarning)

图中给出了：

​

​

因为T大于1e4的值就是问题的开始。在这个值之后，它应该几乎保持不变，但你可以看到它显示的可怕的峰值。我尝试修改epsabs和epsrel，但不起作用。我试着拆分积分限制，显然从E=m到1e4它工作得很好，但当它从E= 1e5积分到无穷大时，又出现了问题。

当我打印积分的结果时，它给出了nn的荒诞值(它可以从5.454320597790705e+17到-3085930898.2363224)。我不知道还能做什么。提前感谢！

Answer 1

您自己提供了答案:整数值变得绝对巨大！这里没有epsabs是有意义的，所以最好通过将它设置为np.inf来禁用它。有了这一点，人们就可以得到合理的

​

​

import matplotlib.pyplot as plt
from scipy.integrate import quad
import numpy as np


m = 100
g = 2.0

TT = np.logspace(-1, 15, 100)
nn = [
    g
    / (2 * np.pi ** 2)
    * quad(
        lambda E: E * np.sqrt(E ** 2 - m ** 2) / (np.exp(E / T) - 1),
        m,
        np.inf,
        epsabs=np.inf,
    )[0]
    for T in TT
]

plt.loglog(TT, nn, "-")
plt.grid()
plt.show()