2

我正在尝试解决以下问题： 在此处输入图像描述

其中，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)

情节给出：

在此处输入图像描述

对于大于 1e4 的 T 值是问题开始的地方。在这个值之后它应该几乎保持不变，但是你可以看到它显示的可怕的峰值。我尝试修改 epsabs 和 epsrel 但没有奏效。我试图拆分积分限制，显然从 E = m 到 1e4 它工作正常，但是当它从 E = 1e5 积分到无穷大时又会出现问题。

当我打印积分的结果时，它给出了 nn 的荒谬值（它可以从 5.454320597790705e+17 变为 -3085930898.2363224）。我不知道还能做什么。提前致谢！

4

1 回答 1

0

您自己已经给出了答案：积分值变得绝对巨大！Noepsabs在这里没有意义，所以最好通过将其设置为来禁用它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()
于 2019-11-13T13:10:42.803 回答