对超过 10^1000 的数字进行无穷大积分

Question

我正在尝试使用包含 M 次幂的项的方程来执行从 - 无穷大到无穷大的积分。在 M 值非常高（5000+）时，该项的值可能超过 10^100 .

我有这个 python (3.10) 程序，它适用于 M = 1000，但除此之外就没有什么了。

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)

我得到的错误是：

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))

integral_value 的结果是

nan

。

有没有一种方法可以使用 Python 执行此积分，这会很方便，因为它允许我将其包含在程序的其余部分中，或者我将被迫使用不同的工具？

Answer 1

假设函数定义正确（考虑到所讨论的幅度，这将是令人惊讶的），不要将其积分到线性空间中。相反，使用对数等价物：

def integrand_log(u, M: float, Z: float, r: float):
    sqrt_r = np.sqrt(r)
    sqrt_term = np.sqrt(2*(1 - r))
    arg1 = (sqrt_r*u + Z)/sqrt_term
    arg2 = (sqrt_r*u - Z)/sqrt_term
    erf_term = special.erf(arg1) - special.erf(arg2)
    return (
        M*np.log(erf_term)
        - 0.5 * (u**2 + np.log(2 * np.pi))
    )

并应用日志身份