A tricky double integral What is
$$\int_0^1 \int_0^1 \frac{ dx \; dy}{1+xy+x^2y^2} ? $$
Can you do one of the integrals and turn it into a single integral?
I get lost in a sea of inverse tangents.
 A: You can rewrite the integral as
$$\int_0^1 \int_0^1 {1 - xy \over 1 - x^3y^3}\,dx\,dy$$
Expanding ${\displaystyle{1 \over 1 - x^3y^3}}$ as a geometric series $\sum_n x^{3n}y^{3n}$
this becomes
$$\int_0^1 \int_0^1 \sum_{n=0}^{\infty} x^{3n}y^{3n} - \sum_{n=0}^{\infty} x^{3n+1}y^{3n+1}\,dx\,dy$$
Integrating this termwise this becomes
$$\sum_{n = 0}^{\infty} {1 \over (3n + 1)^2 } - \sum_{n = 0}^{\infty} {1 \over (3n + 2)^2 } $$
I don't think these sums have a closed form in terms of well-known functions; in effect these are the polygamma functions that Julian Aguirre put in his answer.
A: Integrating first in $y$ gives
$$
-\int_0^1\frac{\pi-6\,\arctan\Bigl(\dfrac{2\,x+1}{\sqrt{3}}\Bigr)}{3\,\sqrt{3} x}\,dx,
$$
which you probably already know. Mathematica does not evaluate it.
The fact that the integrand depends on $x\,y$ suggests the change of variables $\xi=x\,y$, $\eta=y$. The domain $\{(x,y):0\le x\le1,\ 0\le y\le x\}$ is transformed into $\{(\xi,\eta):0\le\xi\le1,\ \xi\le\eta\le\sqrt{\xi}\,\}$, and $dx\,dy=d\xi\,d\eta/\eta$. Then
$$\begin{align*}
\int_{0}^{1} \int_{0}^{1} \frac{ dx\,dy}{1+x\,y+x^2\,y^2}&=2\int_{0}^{1}  \int_{0}^{x} \frac{dx\,dy}{1+x\,y+x^2\,y^2}\\
&=2\int_0^1\int_{\xi}^{\sqrt\xi}\frac{1}{1+\xi+\xi^2}\frac{d\xi\,d\eta}{\eta}\\
&=-\int_0^1\frac{\log\xi}{1+\xi+\xi^2}\,d\xi.
\end{align*}$$
This last integral is not elementary, but at least Mathematica gives a value interms of the PolyGamma function:
$$
\frac{1}{9}\Bigl(\operatorname{PolyGamma}[1, \frac23]-\operatorname{PolyGamma}[1,\frac13]\Bigr).
$$
The numerical value is $0.781302$.
A: The first step would be one $\arctan$ only. 
The denominator is 
$$1+xy+x^2y^2 =\left(\frac{1}{2}+xy\right)^2 + \frac{3}{4}$$
Hence, with $a=\frac{2}{\sqrt{3}}$ we have
$$\int\frac{dy}{1+xy+x^2y^2} =\int\frac{dy}{\left(\frac{1}{2}+xy\right)^2 + 1/a^2}=\qquad\qquad\qquad\text{               }\\
a^2\int\frac{dy}{a^2\left(\frac{1}{2}+xy\right)^2 + 1}=
a^2\frac{\arctan(\frac{a}{2}+axy)}{ax} + C=\\a\frac{\arctan\left(\frac{1}{\sqrt{3}}+\frac{2xy}{\sqrt{3}}\right)}{x} + C$$
so that $$\int_0^1\frac{dy}{1+xy+x^2y^2} = a\frac{\arctan(\frac{1}{\sqrt{3}}+\frac{2x}{\sqrt{3}})}{x} -\frac{\pi}{3\sqrt{3}}\frac{1}{x}$$
The integral of $\arctan(a+bx)/x$ is not that easy, it ends up in dilogarithms... 
