Let us define a generalized mean of two positive real numbers $\mu(a,b)$ as:
$$\frac{1}{2 \mu(a,b)}=\int_0^\infty \frac{dx}{e^{ax}+e^{bx}}=\int_0^1 \frac{dt}{t^{1-a}+t^{1-b}}$$
This integral has a closed form in terms of hypergeometric function or generalized harmonic numbers (note that we can exchange $a$ and $b$ due to symmetry):
$$\frac{1}{2 \mu(a,b)}=\frac{1}{a} {_2 F_1 } \left(1, \frac{a}{a-b};1+\frac{a}{a-b};-1 \right)$$
$$\frac{1}{2 \mu(a,b)}=\frac{1}{2(a-b)} \left( H_{b/(2(a-b))}-H_{(2b-a)/(2(a-b))} \right)$$
Or, using the hypergeometric series, we actually obtain a very simple form for the integral (assuming $a>b$):
$$\frac{1}{2 \mu(a,b)}=\sum_{n=0}^\infty \frac{(-1)^n}{a+(a-b)n}$$
While I'm mainly interested in finding an interated means algorithm to compute the integral in question (like arithmetic-geometric mean for elliptic integrals), I have started with some bounds.
From experiments we have from below:
$$\mu(a,b) \geq \sqrt{\frac{a^2+b^2}{2}} \geq \frac{a+b}{2}$$
And from above:
$$\mu(a,b) \leq \sqrt[3]{\frac{a^3+b^3}{2}} \leq \frac{a^2+b^2}{a+b}$$
The upper bounds are very sharp for $a$ and $b$ close to each other.
Iterating power means for $p=2,3$ we get a mean which intersects $\mu(a,b)$ (red points).
I tried also intermediate power means $\sqrt[p]{\frac{1+x^p}{2}}$ ($p \in (2,3)$), but they mostly intersect $\mu(1,x)$ close to $x=0$.
However, by solving:
$$2 \mu(1,0)=2^{1/p}$$
We find the better bound:
$$\mu(a,b) \geq \sqrt[p]{\frac{a^p+b^p}{2}}$$
Where:
$$p=\frac{\log (2)}{\log (\log (4))}=2.122089644\dots$$
My questions are:
How can we prove that $\mu(a,b) \leq \sqrt[3]{\frac{a^3+b^3}{2}}$?
Can we prove some sharp bound from below as well? Or a sharper bound from above? (In terms of well known means, or at least elementary functions).




This doesn't derive an iterative method to compute the integral, but it does get a different form in terms of the extended Harmonic numbers (although the value is indeed equal).
Preliminary $$ \begin{align} \sum_{k=1}^\infty\frac{(-1)^{k-1}}{k+x} &=\lim_{n\to\infty}\left(\sum_{k=1}^{2n}\frac1{k+x}-2\sum_{k=1}^n\frac1{2k+x}\right)\\ &=\lim_{n\to\infty}\left(\sum_{k=n+1}^{2n}\frac1{k+x}+\sum_{k=1}^n\left(\frac1{k+x}-\frac1{k+x/2}\right)\right)\\ &=\log(2)+\sum_{k=1}^\infty\left(\frac1k-\frac1{k+x/2}\right)-\sum_{k=1}^\infty\left(\frac1k-\frac1{k+x}\right)\\[9pt] &=\log(2)+H_{x/2}-H_x\\[15pt] %&=\log(2)+\psi\!\left(\tfrac x2+1\right)-\psi(x+1) \end{align} $$ where the Extended Harmonic Numbers are related to the digamma function by $H_x=\gamma+\psi(x+1)$
Integral
Assuming $a\gt b$: $$ \begin{align} \int_0^\infty\frac{\mathrm{d}x}{e^{ax}+e^{bx}} &=\int_0^\infty\frac{e^{-bx}\,\mathrm{d}x}{e^{(a-b)x}+1}\\ &=\sum_{k=1}^\infty(-1)^{k-1}\int_0^\infty e^{-(b+k(a-b))x}\mathrm{d}x\\ &=\frac1{a-b}\sum_{k=1}^\infty\frac{(-1)^{k-1}}{\frac{b}{a-b}+k}\\ &=\left\{\begin{array}{} \frac1{a-b}\left(\log(2)+H_{\frac b{2(a-b)}}-H_{\frac b{a-b}}\right)&\text{if $a\gt b$}\\ \frac1{2a}&\text{if $a=b$} \end{array}\right. \end{align} $$
An Identity $$ \begin{align} H_x+H_{x-1/2} &=\sum_{k=1}^\infty\left(\frac1k-\frac1{k+x}\right)+\sum_{k=1}^\infty\left(\frac1k-\frac1{k+x-1/2}\right)\\ &=2\lim_{n\to\infty}\sum_{k=1}^n\left(\frac1k-\frac1{2k+2x}-\frac1{2k-1+2x}\right)\\ &=2\lim_{n\to\infty}\sum_{k=1}^{2n}\left(\frac1k-\frac1{k+2x}\right)-2\lim_{n\to\infty}\sum_{k=n+1}^{2n}\frac1k\\ &=2H_{2x}-2\log(2) \end{align} $$ This shows that my answer and the value in the question are the same.