How to evaluate $\displaystyle\int_{-\infty}^{\infty}e^{-\frac{x^{2}}{2}}\ln\left(1+e^{x}\right)\mathrm{d}x\textbf{ efficiently}$?

Question

I am struggling in evaluating the follow integral: $$\int_{-\infty}^{\infty}e^{-\frac{x^{2}}{2}}\ln\left(1+e^{x}\right)\mathrm{d}x\overset{\text{Wolfram}}{\approx} 2.02049$$ $\color{red}{\text{Despite asking Wolfram to give me 10 digits, it always gives me this result}}$

My process

$$\int_{-\infty}^{\infty}e^{-\frac{x^{2}}{2}}\ln\left(1+e^{x}\right)\mathrm{d}x=1+2\int_{0}^{\infty}e^{-\frac{x^{2}}{2}}\ln\left(1+e^{-x}\right)\mathrm{d}x$$

$$\displaystyle\ln(1+x)=-\sum_{n=1}^{\infty}\frac{(-1)^n}{n}x^n\qquad \text{for }|x|<1$$

$$\begin{align}\int_{-\infty}^{\infty}e^{-\frac{x^{2}}{2}}\ln\left(1+e^{x}\right)\mathrm{d}x=&1-2\sum_{n=1}^{\infty}\frac{\left(-1\right)^{n}}{n}\int_{0}^{\infty}e^{-\frac{x^{2}}{2}}e^{-nx}\mathrm{d}x\\ =&1-2\sum_{n=1}^{\infty}\frac{\left(-1\right)^{n}}{n}\sqrt{\frac{\pi}{2}}e^{\frac{n^{2}}{2}}\text{erfc}\left(\frac{n}{\sqrt{2}}\right)\\ =&1-\sqrt{2\pi}\sum_{n=1}^{\infty}\frac{\left(-1\right)^{n}}{n}e^{\frac{n^{2}}{2}}\text{erfc}\left(\frac{n}{\sqrt{2}}\right)\end{align}$$

This solution is correct, but has an huge defect: it is extremely inefficient to calculate

When I calculate the series for the first $n$ terms I have to take $e^{\frac{n^2}{2}}$ into account

• $n=10$: $e^{\frac{n^2}{2}}\approx 5.18\cdot 10^{21}$ and the series is $\approx \color{blue}{2.0}115$
• $n=20$: $e^{\frac{n^2}{2}}\approx 7.22\cdot 10^{86}$ and the series is $\approx \color{blue}{2.0}181$
• $n=30$: $e^{\frac{n^2}{2}}\approx 2.70\cdot 10^{195}$ and the series is $\approx \color{blue}{2.0}194$
• $n=50$: $e^{\frac{n^2}{2}}\approx \color{red}{7.38\cdot 10^{542}}$ and the series is $\approx \color{blue}{2.020}0$
• $n=110$: $e^{\frac{n^2}{2}}\approx \color{red}{3.03\cdot 10^{2627}}$ and the series is $\approx \color{blue}{2.0204}$

In reality these $n$ are even higher because the series is alternated, with $n+1$ they lose a correct figure and with $n+2$ they recover it.

So to calculate at least $3$ correct decimal digits we need to calculate a number that not even Matlab and Desmos are able to handle (the first Desmos and Matlab overflow number is $2^{1024}\approx 1.79\cdot 10^{308}$)

---

Question

Anyone knows some tricks to solve it in some other way?

Hint 1
The Taylor series of $\ln(1+e^{-x})$ cannot be used, since it has convergence radius $\pi$: $$\ln\left(1+e^{-x}\right)=\ln(2)-\frac{x}{2}-\frac{1}{4}\sum_{n=1}^{\infty}\frac{E_{2n-1}(0)}{n(2n-1)!}x^{2n}\qquad |x|<\pi$$ Therefore used in the integral for $|x|>\pi$ we have that the integral diverges.

• $E_{\nu}(z)$ is the Euler polynomials.

Hint 2
(It's not useful but it's an interesting thing in my opinion)

$$\text{erfc}(z)e^{z^2}\propto \frac{1}{\sqrt{\pi}z}\sum_{n=0}^{\infty}\left(-1\right)^{n}\frac{\left(2n\right)!}{n!}\left(2z\right)^{-2n}$$

Hint 3
I tried even doing the Gaussian series: $$e^{-\frac{x^2}{2}}=\sum_{n=0}^{\infty}\frac{\left(-1\right)^{n}x^{2n}}{2^{n}n!}\qquad \text{for }x\in\mathbb{R}$$ $$\begin{align}\int_{-\infty}^{\infty}e^{-\frac{x^{2}}{2}}\ln\left(1+e^{x}\right)dx=&1+2\int_{0}^{\infty}e^{-\frac{x^{2}}{2}}\ln\left(1+e^{-x}\right)dx\\ =&1+2\int_{0}^{\infty}\sum_{n=0}^{\infty}\frac{\left(-1\right)^{n}x^{2n}}{2^{n}n!}\ln\left(1+e^{-x}\right)dx\\ =&1+2\sum_{n=0}^{\infty}\frac{\left(-1\right)^{n}}{2^{n}n!}\int_{0}^{\infty}x^{2n}\ln\left(1+e^{-x}\right)dx\\ \end{align}$$ $$\color{green}{\int_{0}^{\infty}x^{2n}\ln\left(1+e^{-x}\right)dx=\left(1-2^{-2n-1}\right)\zeta\left(2+2n\right)\left(2n\right)!}$$ $$1+2\sum_{n=0}^{\infty}\frac{\left(-1\right)^{n}}{2^{n}n!}\left(1-2^{-2n-1}\right)\zeta\left(2+2n\right)\left(2n\right)!$$ There must be a mistake here because you end up with a $\color{red}{\text{diverget series}}$

Hint 4
I would like to try the Gauss–Hermite quadrature formula but it seems quite complicated to implement (I guess the result comes, but I don't know how quickly)

Answer 1

There are a number of ways to accelerate the convergence of a given sequence or series. A common tool in the context of alternating series is the Euler transform, which states that $$ \sum_{n=0}^{\infty} (-1)^n a_n = \sum_{n=0}^{\infty} (-1)^n \frac{\Delta^n a}{2^{n+1}}, $$ where $$ \Delta^n a = \sum_{k=0}^{n} (-1)^n \binom{n}{k} a_{n-k}. $$

To apply that to your approach, we would take $$ a_n = \frac{1}{n}e^{\frac{n^{2}}{2}}\text{erfc}\left(\frac{n}{\sqrt{2}}\right). $$

Using SageMath we can implement this as follows:

n,k = var('n,k')
a(n) = exp(n^2/2)*erfc(n/sqrt(2))/n
def fd(n):
  return sum([(-1)^k * binomial(n,k)*a(n-k) for k in range(n)])
sum_result = N(1 - sqrt(2*pi) * sum([(-1)^n * fd(n)/2^(n+1) for n in range(1,50)]))

We can also work out the integral numerically:

int_result = numerical_integral(exp(-x^2/2)*log(1+exp(x)), -10,10)

Given that $\log(1+e^x) \sim x$ for large $x$, it's pretty easy to show should be within machine precision of the improper integral.

You can see how well this worked in this Sage Cell. There we see that

sum_result = 2.02049074000466
int_result = (2.020490740004666, 2.2431953213344902e-14)

Note that the integral result includes an error estimate and the sum result is quite close.

Of course, the numerical check rather begs the question - if numerical computation of the integral is your main objective, why not just use numerical integration techniques in the first place?