Четырехзначный интеграл с бесконечным пределом - PullRequest
Новая
       88

Четырехзначный интеграл с бесконечным пределом

2 голосов
/ 08 марта 2019

Я пытаюсь решить следующее: enter image description here

где 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)

А сюжет дает:

enter image description here

Для значений T больше 1e4 начинается проблема. После этого значения оно должно оставаться почти постоянным, но вы можете видеть ужасные пики, которые он показывает. Я попытался изменить epsabs и epsrel, но это не сработало. Я пытался разделить пределы интегрирования, и, по-видимому, от E = m до 1e4 он работает нормально, но когда он интегрируется от E = 1e5 до бесконечности, то снова возникают проблемы.

Когда я печатаю результаты интеграла, он дает абсурдные значения для nn (он может идти от 5.454320597790705e + 17 до -3085930898.2363224). Я не знаю, что еще делать. Заранее спасибо!

...