A possibility to use the fact that $F$ is homogeneous is to pass to polar coordinates, $x=r\; \cos t$, $y=r\; \sin t$, so
$$dx\wedge dy = (\cos t\; dr - r\; \sin t\; dt)\wedge (\sin t\; dr + r\; \cos t\; dt) =r\; dr\wedge dt\ .$$
If $F$ is homogeneous of degree $k$, then we have:
$$ \int_{B_1(0)}F(x,y)\; dx\; dy = \int_0^1\int_0^{2\pi} r^k\; F(\cos t, \sin t)\; r\; dr\; dt = \frac 1{k+2}\int_0^{2\pi} F(\cos t, \sin t)\; dt\ .$$
For example, for the function $F(x,y)=x^8+y^8$ we have $k=8$ and we can compute exactly, and numerically:
var( 'x,y,t' )
F(x,y) = x^8 + y^8
k = 8
J = 1/(k+2) * integrate( F( cos(t), sin(t) ), (t,0,2*pi) )
print "J = %s ~ %s" % ( J, J.n() )
This gives:
J = 7/64*pi ~ 0.343611696486384
Using the idea of nbruin:
sage: integral( integral( F(x,y), (y,-sqrt(1-x^2),+sqrt(1-x^2)) ), (x,-1,1) )
7/64*pi
The transformation to polar coordinates also works for a "positively homogeneous" function like in the following example:
F(x,y) = abs(x)^(3/5) * abs(y)^(7/5)
k = 2
J = 1/(k+2) * integrate( F( cos(t), sin(t) ), (t,0,2*pi) )
print "J is %s" % J
print "J is numerically %s" % J.n()
Result:
J is 1/4*integrate(abs(cos(t))^(3/5)*abs(sin(t))^(7/5), t, 0, 2*pi)
J is numerically 0.534479668623671
The integral could not be computed, but the numerical value is accessible.
Alternatively, we can use the numerical integral, the following example shows which is the error, the trap i am falling in every time:
sage: x=0.12345; numerical_integral( lambda y: F(x,y), (-sqrt(1-x^2),+sqrt(1-x^2)) )
(0.23319001033141668, 1.3392247519283806e-09)
(My error is that i always forget about the error... The result of the numerical_integral
is a tuple, not a number.) In order to get a number for this particular value of $x$, we have to take the pythonically zeroth part of the tuple.
sage: x=0.12345; numerical_integral( lambda y: F(x,y), (-sqrt(1-x^2),+sqrt(1-x^2)) )[0]
0.23319001033141668
So the numerical integral using Fubini is:
sage: numerical_integral( lambda x: numerical_integral( lambda y: F(x,y), (-sqrt(1-x^2),+sqrt(1-x^2)) )[0], (-1,1) )[0]
0.5344796700960925
Doesn't integral(integral(F,x,-sqrt(1-y^2),sqrt(1-y^2)),y,-1,1) do the trick? You shouldn't expect too nice an answer from this, of course: for F=1 the value of the integral is already a transcendental number. The standard calculus bag of tricks applies, so depending on F there may be smart substitutions that might get you a closed-form solution that is otherwise not easily attained. Rewriting in polar coordinates would be an obvious thing to try too.