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

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)

• change variable with the logarithm Apr 26, 2023 at 22:37
• Hint. Split the integral over the domain in two: over the negative domain and over the positive one. In the negative domain, you do have that $e^x<1.$ For the positive domain, factor out $e^x$ from the logarithm, to make -- in your words -- the argument of the logarithm "finite". Apr 26, 2023 at 22:47
• Rewrite the logarithm’s argument as $e^x\cdot (1 + e^{-x})$, then expand $\log(1 + e^{-x})$ in powers of $e^{-x}$ for the $x>0$ case, etc. Apr 26, 2023 at 22:57
• Mathematica can go to more digits of precision using NIntegrate and WorkingPrecision, e.g., NIntegrate[Exp[-x^2/2] Log[1 + E^x], {x, -Infinity,Infinity}, WorkingPrecision -> 20] gives 2.0204907400046658800. Jan 5 at 20:59
• How about using dlmf.nist.gov/7.9 ?
– Gary
Jan 11 at 23:21

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?