Как решить проблему ошибки округления в scipy.integrate.quad для конкретного интегранда? ⇐ Python

[Гость]

Я провожу следующий численный расчет.

Код: Выделить всё

import numpy as np
from scipy.integrate import quad

def E_k(gap, q, mu):
return np.sqrt(abs(gap) * abs(gap) + (q * q - mu) * (q * q - mu))

quad(lambda q: q * q * (1 - (q * q - mu) / E_k(gap, q, mu) * np.tanh(E_k(gap, q, mu) / 2 / T)),
0, np.inf, epsabs=1.49e-6, epsrel=1.49e-6, limit=50)

It is aimed to set the value of mu to be in the range of (-30, 1), the value of gap to be in the range of (0, 3), and T to be positive. The quad will give wrong results (it tells me roundoff error detected) when T is small (approximately less than 0.1) and mu is positive or negatively close to zero (approximately -0.0005 < mu < 1).
I have tried to divide the integrating interval at point mu when mu is positive.

Код: Выделить всё

quad(lambda q: q * q * (1 - (q * q - mu) / E_k(gap, q, mu) * np.tanh(E_k(gap, q, mu) / 2 / T)),
0, mu, epsabs=1.49e-6, epsrel=1.49e-6, limit=50)[0] + \
quad(lambda q: q * q * (1 - (q * q - mu) / E_k(gap, q, mu) * np.tanh(E_k(gap, q, mu) / 2 / T)),
mu, np.inf, epsabs=1.49e-6, epsrel=1.49e-6, limit=50)[0]

But it seems not to work. I have also increased the limit and adjusted epsabs, epsrel in quad, but the results are still uncorrected.


Источник: https://stackoverflow.com/questions/759 ... cific-inte