Is there an elegant way to evaluate $ I={ \int \sqrt[8]{\frac{x+1}{x}} \ \mathrm{d}x}$?

Just to rework @H. R.'s solution by hand. First off setting $$t=\frac{x+1}x\ge0$$ was definitely a good step. Solving for $x$, $$x=\frac1{t-1}$$ Note that this implies that either $0\le t<1$ or $t>1$, a property which we will use later to simplify the result. Then $$\int\sqrt[8]{\frac{x+1}x}dx=-\int\frac{t^{1/8}}{(t-1)^2}dt=-8\int\frac{u^8}{(u^8-1)^2}du$$ Where we have let $t=u^8$, just as @H. R. Further, observe that
$$\frac d{du}\left(\frac u{u^8-1}\right)=\frac{-7u^8-1}{(u^8-1)^2}$$ So that $$\int\sqrt[8]{\frac{x+1}x}dx=\int\left[\frac d{du}\left(\frac u{u^8-1}\right)-\frac1{u^8-1}\right]du=\frac u{u^8-1}-\int\frac{du}{u^8-1}=\frac u{u^8-1}-J$$ Then partial fractions is easy: $$\frac1{u^8-1}=\frac1{\prod_{k=0}^7(u-\omega_k)}=\sum_{k=0}^7\frac{A_k}{u-\omega_k}$$ Where $\omega_k=e^{i\theta_k}=\cos\theta_k+i\sin\theta_k$ and $\theta_k=\frac{\pi k}4$. Evaluating $$\lim_{u\rightarrow\omega_k}\frac{u-\omega_k}{u^8-1}=\lim_{u\rightarrow\omega_k}\frac1{8u^7}=\frac1{8\omega_k^7}=\frac{\omega_k}8=\lim_{u\rightarrow\omega_k}\sum_{j=0}^7A_j\frac{u-\omega_k}{u-\omega_j}=\sum_{j=0}^7A_j\delta_{jk}=A_k$$ So $$\begin{align}J&=\frac18\sum_{k=0}^7\int\frac{\omega_k}{u-\omega_k}du=\frac18\int\left[\frac1{u-1}-\frac1{u+1}+\sum_{k=1}^3\frac{2u\cos\theta_k-2}{u^2-2u\cos\theta_k+1}\right]du\\ &=\frac18\int\left[\frac1{u-1}-\frac1{u+1}+2\sum_{k=1}^3\frac{(u-\cos\theta_k)\cos\theta_k-\sin^2\theta_k}{(u-\cos\theta_k)^2+\sin^2\theta_k}\right]du\\ &=\frac18\left[\ln\left|\frac{u-1}{u+1}\right|+\sum_{k=1}^3\left(\cos\theta_k\ln\left(u^2-2u\cos\theta_k+1\right)-2\sin\theta_k\tan^{-1}\left(\frac{u-\cos\theta_k}{\sin\theta_k}\right)\right)\right]+C\end{align}$$ We may apply the values $\cos\theta_1=-\cos\theta_3=\sin\theta_1=\sin\theta_3=\frac1{\sqrt2}$, $\cos\theta_2=0$, and $\sin\theta_2=0$ to arrive at $$\begin{align}J&=\frac18\left[\ln\left|\frac{u-1}{u+1}\right|+\frac1{\sqrt2}\ln\left(\frac{u^2-\sqrt2u+1}{u^2+\sqrt2u+1}\right)\right.\\ &\left.-\sqrt2\tan^{-1}\left(\sqrt2u-1\right)-\sqrt2\tan^{-1}\left(\sqrt2u+1\right)-2\tan^{-1}u\right]+C\end{align}$$ Now we can add two of those arctangents because $u$ can't take on values that make the new denominator $0$ nor change the sign of the numerator when the denominator is negative: $$\begin{align}J&=\frac18\left[\ln\left|\frac{u-1}{u+1}\right|+\frac1{\sqrt2}\ln\left(\frac{u^2-\sqrt2u+1}{u^2+\sqrt2u+1}\right)-\sqrt2\tan^{-1}\left(\frac{\sqrt2u}{1-u^2}\right)-2\tan^{-1}u\right]+C\end{align}$$ We now have the solution

$$\begin{align}\int\sqrt[8]{\frac{x+1}x}dx&=\frac{\sqrt[8]{\frac{x+1}x}}{\frac{x+1}x-1}-\frac18\left[\ln\left|\frac{\sqrt[8]{\frac{x+1}x}-1}{\sqrt[8]{\frac{x+1}x}+1}\right|+\frac1{\sqrt2}\ln\left(\frac{\sqrt[4]{\frac{x+1}x}-\sqrt2\sqrt[8]{\frac{x+1}x}+1}{\sqrt[4]{\frac{x+1}x}+\sqrt2\sqrt[8]{\frac{x+1}x}+1}\right)\right.\\ &\left.-\sqrt2\tan^{-1}\left(\frac{\sqrt2\sqrt[8]{\frac{x+1}x}}{1-\sqrt[4]{\frac{x+1}x}}\right)-2\tan^{-1}\sqrt[8]{\frac{x+1}x}\right]+C\end{align}$$

Which matches the Mathematica solution.

I have Mathematica 10.4 and what it returns for the primitive function is too lengthy! Off the top of my head, I don't think that there is any other way but the strategy to obtain the final answer is not that much complicated.! :) I just explain the path that you may go to obtain this result.

Firstly, we set

$$\begin{align} t &= \frac{x+1}{x} \\ dx &= \frac{1}{(t-1)^2} dt \end{align}$$

and the integral will become

$$I=-\int \frac{\sqrt[8]{t}}{(t-1)^2} dt $$

and then setting

$$\begin{align} t &= u^8 \\ dt &= 8u^7du \end{align}$$

will lead to

$$I=-8\int \frac{u^8}{(u^8-1)^2} du$$

and from here on, you should use partial fraction to reach what Mathematica gives!

Through the substitution $1+\frac{1}{x}=u$, followed by $u=v^8$, we get:

$$ \int \sqrt[8]{1+\frac{1}{x}}\,dx = \int\frac{u^{1/8}}{(1-u)^2}\,du =\int\frac{8v^8\,dv}{(1-v^8)^2}$$ and the last integral can be computed through partial fraction decomposition (obviously, the eighth roots of unity are involved).