I am attempting to perform on an integral from -infinity to infinity, with an equation containing a term that is raised to the power M. At very high (5000+) values of M, this value of this term can exceed 10^100.
I have this python (3.10) program, which works for M = 1000, but nothing much beyond that.
import numpy as np
from scipy import integrate, special
def integrand(u, M, Z, r):
sqrt_r = np.sqrt(r)
sqrt_term = np.sqrt(2 * (1 - r))
arg1 = (Z + sqrt_r * u) / sqrt_term
arg2 = (-Z + sqrt_r * u) / sqrt_term
erf_term = special.erf(arg1) - special.erf(arg2)
return np.exp(-u**2 / 2) / np.sqrt(2 * np.pi) * erf_term**M
def evaluate_integral(M, Z, r):
result, _ = integrate.quad(integrand, -np.inf, np.inf, args=(M, Z, r))
return 1 / (2**M) * result
# Example usage:
M = 1000
Z = 5
r = 0.6
integral_value = evaluate_integral(M, Z, r)
The error I get is:
IntegrationWarning: The occurrence of roundoff error is detected, which prevents
the requested tolerance from being achieved. The error may be
underestimated.
result, _ = integrate.quad(integrand, -np.inf, np.inf, args=(M, Z, r))
and the result for integral_value is nan.
Is there a way I can perform this integral using Python, which would be convenient as it allows me to be included with the rest of my program, or will I be forced to use a different tool?
No, your only problem is that you took the factor 1/2**M outside of the integral. If you put it back inside you will then have
(erf_term/2)**M
and that will not blow up because the difference of two error functions cannot exceed 2, so the new bracketed factor is less than 1.