Answer
$\mathop \smallint \limits_0^1 \mathop \smallint \limits_0^1 y\sqrt {1 + xy} {\rm{d}}y{\rm{d}}x \simeq 0.575161$
Work Step by Step
We use Fubini's Theorem to change the order of integration and re-write the integral:
$I = \mathop \smallint \limits_{x = 0}^1 \mathop \smallint \limits_{y = 0}^1 y\sqrt {1 + xy} {\rm{d}}y{\rm{d}}x = \mathop \smallint \limits_{y = 0}^1 \left( {\mathop \smallint \limits_{x = 0}^1 y\sqrt {1 + xy} {\rm{d}}x} \right){\rm{d}}y$
We evaluate first the inner integral with respect to $x$.
$S\left( y \right) = \mathop \smallint \limits_{x = 0}^1 y\sqrt {1 + xy} {\rm{d}}x$
The antiderivative of $y\sqrt {1 + xy} $ with respect to $x$ is $\frac{2}{3}{\left( {1 + xy} \right)^{3/2}}$ because
$\frac{\partial }{{\partial x}}\left( {\frac{2}{3}{{\left( {1 + xy} \right)}^{3/2}}} \right) = y\sqrt {1 + xy} $
So,
$S\left( y \right) = \mathop \smallint \limits_{x = 0}^1 y\sqrt {1 + xy} {\rm{d}}x = \frac{2}{3}{\left( {1 + xy} \right)^{3/2}}|_{x = 0}^1$
$ = \frac{2}{3}{\left( {1 + y} \right)^{3/2}} - \frac{2}{3}$
Thus,
$I = \mathop \smallint \limits_{y = 0}^1 \left( {\mathop \smallint \limits_{x = 0}^1 y\sqrt {1 + xy} {\rm{d}}x} \right){\rm{d}}y$
$ = \mathop \smallint \limits_{y = 0}^1 \left( {\frac{2}{3}{{\left( {1 + y} \right)}^{3/2}} - \frac{2}{3}} \right){\rm{d}}y$
$ = \frac{2}{3}\mathop \smallint \limits_{y = 0}^1 {\left( {1 + y} \right)^{3/2}}{\rm{d}}y - \frac{2}{3}\mathop \smallint \limits_{y = 0}^1 {\rm{d}}y$
$ = \frac{4}{{15}}{\left( {1 + y} \right)^{5/2}}|_0^1 - \frac{2}{3}$
$ = \frac{{16}}{{15}}\sqrt 2 - \frac{4}{{15}} - \frac{2}{3}$
$ = \frac{{16}}{{15}}\sqrt 2 - \frac{{14}}{{15}} \simeq 0.575161$
So, $\mathop \smallint \limits_0^1 \mathop \smallint \limits_0^1 y\sqrt {1 + xy} {\rm{d}}y{\rm{d}}x \simeq 0.575161$.
