[SOLVED] Performing integral over infinity with numbers in excess of 10^1000

Performing integral over infinity with numbers in excess of 10^1000

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?

Solution

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.