numerical integration - NIntegrating within an Ellipsoid

I need to numerically integrate an expensive positive-definite function over a 2D domain. I know by other ways that the function is basically zero for values outside the following ellipse:

Ellipsoid[center,Sqrt[Eigenvalues[matrix]], Eigenvectors[matrix]]

How do I tell Mathematica to integrate only within the ellipsoid?

---

Here is what I got using your nice suggestions (timings are normalized to the slowest):

100 sec - integration over the rectangle circumscribing the ellipse

55 sec - using Boole to integrate only over ellipse (@b.gatessucks)

47 sec - directly limiting the domain to the ellipse (@whuber)

41 sec - mapping the ellipse to a unit circle (@Jens)

Answer

Why not use elliptical coordinates? That's what they are there for.

For example, if your function is f[x_,y_], then you define the coordinate transformation

x[u_, v_]:= Cos[v] Cosh[u];
y[u_, v_]:= Sin[v] Sinh[u];

The Jacobian is Sin[v]^2 + Sinh[u]^2, so you simply do the integral like this:

NIntegrate[
 f[x[u, v], y[u, v]] (Sin[v]^2 + Sinh[u]^2), {u, 0, 1}, {v, 0, 2 Pi}]

That's the basic idea. You can always rotate and scale your function arguments to fit within this particular ellipse which has axis lengths of Cosh[1] and Sinh[1], respectively. You can also adjust the integration limit in {u, 0, 1} to something other than 1 depending on the eccentricity of the ellipse. These steps are easy once you understand the coordinate system.

You can also do all this using the VectorAnalysis package, but I'd say that is overkill for this relatively simple problem.

Rescaling to a circle in 2D

In fact, you may also want to look into just rescaling your x and y coordinates to fit into a circle and use polar coordinates... (that also generalizes more easily to higher dimensions).

There isn't exactly one built-in function to do the rescaling, but you can probably find all you need by starting at this post. I'll just describe the logic:

Assume you have already figured out the semimajor and -minor axes (vec1 and vec2) as well as the center (in a variable center) of the ellipse. Then you can figure out the rotation angle required to make vec1 the x-axis, leading to a rotation R. The scale transformation T that you need is given by the lengths of vec1 and vec2, resp.

Now compose your function f[x,y] with these transformations. This transformation introduces a Jacobian given by the product Norm[vec] * Norm[vec2].

The final step is to go to polar coordinates which introduces an additional factor r (the radial coordinate).

Polar coordinates are probably best in this case, unless you plan on doing other manipulations with your function that go beyond integration.

Edit

Turns out it's almost more difficult to describe the transformations in words than it is to do it in Mathematica. So I'm adding the code now:

rescaling = Composition[TranslationTransform[center], 
   RotationTransform[{{1, 0}, vec1}], 
   ScalingTransform[{Norm[vec1], Norm[vec2]}]]

This is now used in the integral of your function which I assume to be defined as

f[{x_, y_}] := ...

(note that the argument is a single list representing the point). With this, it only remains to write the integral

integrand = Composition[f, rescaling];
polarIntegrand[r_?NumericQ, a_?NumericQ] := 
     integrand[{r Cos[a], r Sin[a]}] r;
NIntegrate[
  polarIntegrand[r,a], 
   {r, 0, 1}, {a, 0, 2 Pi}] Norm[vec1] Norm[vec2]

Here, I've tweaked the integrand a little to make sure we don't get slowed down by Mathematica attempting to symbolically simplify it. The ?NumericQ rules that out. That should be all you need.