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.