% Variable integral Gaussian $\int_{-\infty}^\infty \frac{e^{-x^2}}{1+a e^{-x}} dx$

Question

He estado tratando de calcular esta integral $$\int_{-\infty}^{\infty} \frac{e^{-x^2}}{1+a e^{-x}} dx$$ rápidamente y con un alto grado de exactitud.

Tengo algunos resultados parciales, por ejemplo, para $n \in \mathbb N$ $$ \frac{1}{\sqrt\pi} \int_{-\infty}^{\infty} \frac{e^{-x^2}}{1+e^{-(x+n/2)}}dx = \frac{1}{2}\sum_{i=0}^{2n}(-1)^{i} e^{-i(2n-i)/4},$$ y para $a \in (0,1)$ hemos $$ \frac{1}{\sqrt\pi} \int_{0}^\infty \frac{e^{-x^2}}{1+ae^{-x}}dx = \frac{1}{2}\sum_{n=0}^\infty (-1)^n a^n e^{n^2/4} \mathrm{erfc}(n/2)$$ donde $\mathrm{erfc}$ es el error de la función de complementar el fuction $$ \mathrm{erfc}(x) = 1-\mathrm{erf}(x),$$ la que se obtiene mediante la sustitución de $$ \frac{1}{1+ae^{-x}} = \sum_{n=0}^\infty (-1)^n a^n e^{-nx}.$$ Esta $\mathrm{erfc}$ serie no converge muy rápidamente, a menos que $a$ es muy pequeña, y wolframalpha puede calcular la integral de forma muy precisa entre decir $-100$ $100$ para los distintos valores de $a$, por lo que me debe de faltar un truco.

Editar He encontrado una serie infinita para el conjunto de la integral, pero converge lentamente $$\frac{1}{\sqrt\pi}\int_{-\infty}^\infty \frac{e^{-x^2}}{1+e^{-(x+a)}}dx = \frac{e^{-a^2}}{2}\left[\sum_{n=0}^\infty (-1)^n \mathrm{erfcx}(-a+n/2) + \sum_{n=1}^\infty (-1)^{n+1} \mathrm{erfcx}(a+n/2)\right],$$ donde $\mathrm{erfcx}$ es de la escala de complemento función de error $$ \mathrm{erfcx}(x) = e^{x^2} \mathrm{erfcx}(x) = e^{x^2}(1-\mathrm{erf}(x)).$$

Answer 1

Desde su integrando es en forma de campana y no hay singularidades de $a>0$, luego de Euler-Maclaurin de expansión de la muestra que el simple trapezoidal de cuadratura converge exponencialmente rápido.

El uso de este código, que fue capaz de evaluar su integral a flotante de precisión en 17 integrando las evaluaciones, en el doble en 33 evaluaciones, de largo el doble de 65 años, y de 50 dígitos binarios en 260. Se ha duplicado el número correcto de dígitos por reducir a la mitad el $h$ que se observa.

Suponiendo que el integrando se lleva a 100 ns a evaluar, esto implica que usted puede calcular su integral para doble precisión en ~3.3 microsegundos.

Aquí está el código que he usado

Real a = 2;
auto f = [&](Real t) { return exp(-t*t)/(1+a*exp(-t)); };
Real t_max = sqrt(-log(std::numeric_limits<Real>::epsilon()));
Real Q = boost::math::adaptive_trapezoidal(f, -t_max, t_max, tol);
std::cout << "Q = " << Q << std::endl;

Que presumiblemente está utilizando el entorno de programación, pero la regla trapezoidal es lo suficientemente simple para el código que esto no debería ser un problema.