Proof of a Ramanujan Integral

You can use residues, let $f(z)=\frac{1}{(1+z^2)(1+r^2z^2)(1+r^4z^2)...}$ this has an infite set of singularities at $z=\pm i 1/r^{n},n\in\Bbb N_0$. We can see that $\large\int_0^\infty\frac{1}{(1+x^2)(1+r^2x^2)(1+r^4x^2)...}dx=\frac{1}{2}\int_{-\infty}^\infty\frac{1}{(1+x^2)(1+r^2x^2)(1+r^4x^2)...}dx$

So we will consider the upper semi circle contour, spanning over $-\infty$ to $+\infty$, only the positive singularities are in this contour, i.e. $z=i\frac{1}{r}$.

Now $\operatorname{Res}(f(z),z=\frac{i}{r^{n}})=\large\lim_{z\to i\frac{1}{r^{n}}}(z-i\frac{1}{r^{n}})f(z)=\large\lim_{z\to i\frac{1}{r^{n}}}\frac{1}{r^n}(r^nz-i)f(z)$

$=\lim_{z\to i\frac{1}{r^{n}}}\frac{1}{r^n(r^nz+i)}\prod_{j=0,j\ne n}^\infty\frac{1}{(1+x^{2}r^{2j})}$

$=\large\frac{1}{2ir^n}\prod_{j=0,j\ne n}^\infty\frac{1}{(1-r^{2j-2n})}$

Now $\frac{1}{2}\int_{-\infty}^\infty\frac{1}{(1+x^2)(1+r^2x^2)(1+r^4x^2)...}dx=\frac{1}{2}2\pi i\sum \operatorname{Res}(f)=\frac{1}{2}\pi i\sum_{n=0}^\infty\frac{1}{ir^n}\prod_{j=0,j\ne n}^\infty\frac{1}{(1-r^{2j-2n})}$

From here I'm not sure how to reach the closed form, but hopefully this helps to show a different approach, even if you are not too well versed in integrals by residue :)

I find this relationship charming.

Here is what I have so far as an alternative solution without resorting to complex analysis. For $n=0,1,...$ define \begin{align} I_n=\int_0^{\infty} \frac{dx}{(1+x^2)(1+r^2x^2)...(1+r^{2n}x^2)} \end{align} Doing many partial fractions I could establish for $n=1,2,3,4,...$ that \begin{align} I_{n}&=\frac{1-r^{2n-1}}{1-r^{2n}}I_{n-1}\\ I_0&=\frac{\pi}{2} \end{align} Therefore \begin{align} I_{\infty}&=\frac{\prod_{n=1}^{\infty}\Big(1-r^{2n-1}\Big)}{\prod_{n=1}^{\infty}\Big(1-r^{2n}\Big)}\frac{\pi}{2} \end{align} Now using the same arguments as in this post, we can establish $(1)$.

Here is a way you can do it with brute force with the aid of change of variables and Fubini's Theorem. Let \begin{align*} I &= \int_{0}^\infty \frac{1}{(1+x^2)(1+r^2x^2) \ \dots (1+r^{2n}x^2)} \ dx. \end{align*} Note for each $i \in \lbrace 1, \ \dots \ ,n \rbrace,$ we can write \begin{align*} \frac{1}{1+r^{2i} x^2} &= \int_{0}^{\infty} e^{-(1+r^{2i}x^2) y_i} \ dy_i. \end{align*} Thus, using Fubini's Theorem,

\begin{align*} I &=\int_{0}^\infty \int_{0}^\infty \ \dots \ \int_{0}^\infty e^{-y_1(1+x^2)} e^{-y_2(1+r^2x^2)} \ \dots \ e^{-y_{n}(1+r^{2n}x^2)} \ dy_n \ \dots \ dy_1\ dx \\ &=\int_{0}^\infty \int_{0}^\infty \ \dots \ \int_{0}^\infty e^{-(y_1 + \ \dots \ + y_n)} \ e^{-(y_1+r^2y_2 \ \dots \ + r^{2n}y_n)x^2} \ dy_n \ \dots \ dy_1 \ dx \\ &=\int_{0}^\infty \int_{0}^\infty \ \dots \ \int_{0}^\infty e^{-(y_1 + \ \dots \ + y_n)} \ e^{-(y_1+r^2y_2 \ \dots \ + r^{2n}y_n)x^2} \ dx \ dy_n \ \dots \ dy_1. \end{align*} To carry out the integral with respect to $x,$ use the one dimensional change of variables $x=\frac{t}{\sqrt{y_1+r^2y_2 \ \dots \ + r^{2n}y_n}}$ and the well-known Gaussian integral $\int_{0}^\infty e^{-t^2} \ dt= \frac{\sqrt{\pi}}{2}.$ We get

\begin{align*} I &= \frac{\sqrt{\pi}}{2} \int_{0}^{\infty} \ \dots \ \int_{0}^{\infty} \frac{e^{-(y_1 + \ \dots \ + y_n)}}{\sqrt{y_1+r^2y_2 \ \dots \ + r^{2n}y_n}} \ dy_n \ \dots \ dy_1. \end{align*}

Now, perform the multidimensional change of variables \begin{align*} y_i &= t_i - r^2 t_{i+1}, \quad 1 \leq i <n \\ y_n &= t_n. \end{align*} By induction on $n,$ it is easy to verify $\left|\frac{\partial(y_1 ,\ \dots \ , y_n)}{\partial(t_1 , \ \dots \ , t_n)}\right|=1$ and that transformed region of integration is

\begin{align*} t_1 &>0 \\ 0 &< t_i < \frac{t_{i-1}}{r^2},\quad 2 \leq i \leq n. \end{align*}

With this, \begin{align*} I &= \frac{\sqrt{\pi}}{2} \int_{0}^{\infty} \int_{0}^{t_1/r^2} \ \dots \ \int_{0}^{t_{n-1}/r^2} \frac{e^{-t_1-(1-r^2)(t_2 + \ \dots \ + t_n)}}{\sqrt{t_1}}\ dt_n \ \dots \ dt_2 \ dt_1. \end{align*}

This makes $I$ an "elementary" integral that can be evaluated either by hand or with Mathematica. To carry out the integration with respect to $t_n, \ \dots \ , t_2,$ one must repeatedly apply the exponential identity $$\int_{0}^{t} e^{-ax} \ dx = \frac{1-e^{-at}}{a}$$ for $a>0,$ which is pretty messy. To carry out the final integration with respect to $t_1,$ make the substitution $t_1=u_1^2$ which will transform the resulting integrand into a Gaussian function.
At the moment, I do not know of a clean, slick way to evaluate this multidimensional integral and arrive at the general answer, but I will look into it.