Numerical integration: $\int_0^1 \sqrt{\int_0^1 f(x,y)dx}\ dy$ Recently, I encountered the following integration:
$$\int_0^1 \sqrt{\int_0^1 f(x,y)dx}\ dy$$
where $f(x,y)$ is rather complex. We know that with the help of Mathematica it is easy to compute the numerical value of $\int_0^1 \int_0^1 f(x,y)dxdy$. However, it becomes useless to handle with the integration I mentioned above. I am wondering if there are any methods that can simplify this integration to the form of $\int_0^1 \int_0^1 g(x,y)dxdy$, or can we find any efficient way to compute its numerical value?
 A: You could define the function $g(y) = \sqrt{\int_0^1 f(x,y)\,dx}$ and then use a standard quadrature method to compute $\int_0^1 g(y)\,dy$. Each time the quadrature routine evaluates $g$ it has to call a separate quadrature routine to evaluate $\int_0^1 f(x,y)\,dx$ with $y$ fixed. I realize this may not be the most efficient method, but it should work. Perhaps with more information about $f(x,y)$ more could be said.
A: As mentioned in the first answer, you first define a function 
$$F(y) = \sqrt{\int_0^1 f(x,y)\,dx}$$ (I guess this will be a numerical procedure, really) and then compute 
$$\int_0^1 F(y) dy.$$ 
From there, I guess it's really a programming question.  To do this in Mathematica (which you mention you have), you would do something like so:
f[x_, y_] := x^2 + y^2;
F[y_?NumericQ] := Sqrt[NIntegrate[f[x, y], {x, 0, 1}]]
NIntegrate[F[y], {y, 0, 1}]

(* Out: 0.796843 *)

I assume this example is much simpler than yours but, hopefully, illustrates the idea.  Also, this can be evaluated in closed form so we can compare:
Integrate[Sqrt[Integrate[f[x, y], {x, 0, 1}]], {y, 0, 1}]

(* Out: 1/Sqrt[3] + ArcSinh[Sqrt[3]]/6 *)

