13
$\begingroup$

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)

Improve this question
$\endgroup$
3
  • $\begingroup$ possible duplicate of Conditional numerical integration boundaries $\endgroup$
    – rm -rf ♦
    Mar 12, 2012 at 17:28
  • 1
    $\begingroup$ @R.M The questions are awfully similar, but IMO there is a subtle difference. The question here is about restricting calculations to a certain area whereas the other question is about totally preventing evaluation whenever the higher bounds are smaller than the lower bounds. $\endgroup$
    – Sjoerd C. de Vries
    Mar 12, 2012 at 18:03
  • $\begingroup$ Could you time, too, for comparison, Jens's method of using elliptical coordinates? Having a complete comparison of all the methods would be edifying. $\endgroup$
    – murray
    Mar 13, 2012 at 19:35

3 Answers 3

Reset to default
11
$\begingroup$

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.

Improve this answer
$\endgroup$
10
  • 1
    $\begingroup$ NIntegrate also does some sort of coordinate transformation to match the shape of the domain specified with Boole (check the points where it evaluates the integral using the code from Heike's comment). Also, it seems to be smart enough to detect any symmetries of the integrand and use them to reduce the number of points. $\endgroup$
    – Szabolcs
    Mar 13, 2012 at 16:56
  • 1
    $\begingroup$ Well, there isn't always a suitable coordinate transformation, so Boole is a good general approach. But elliptical domains are special, and one should exploit that whenever possible. Not knowing how clever Boole is, I instead trust my own cleverness... $\endgroup$
    – Jens
    Mar 13, 2012 at 17:04
  • $\begingroup$ do you think it is gonna be faster in this way? $\endgroup$
    – Valerio
    Mar 13, 2012 at 18:55
  • $\begingroup$ is there a built-in way to change coordinates in the case of rotated and translated ellipsoid? I like the idea of mapping the ellipsoid to a sphere $\endgroup$
    – Valerio
    Mar 13, 2012 at 18:59
  • $\begingroup$ I'll try to add some more to this post to answer this question. $\endgroup$
    – Jens
    Mar 13, 2012 at 19:52
12
$\begingroup$

If you know the equation defining your ellipsoid you could use Boole[] to constrain the integration domain :

myF[x_,y_]=Abs[x+y]
NIntegrate[Boole[(x/3)^2 + (y/2)^2 <= 1] myF[x,y], {x, -5, 5}, {y, -5, 5}]

Note that this will actually prevent myF[x, y] from being evaluated outside the domain specified by Boole. This feature of NIntegrate is described here.

Improve this answer
$\endgroup$
7
  • $\begingroup$ This really doesn't help reducing calculation time as asked by the OP. In fact, adding Boole increases calculation time 6-fold. This is logical as adding Boole does not prevent the evaluation of the expensive function outside of the area where the argument of Boole is true. $\endgroup$
    – Sjoerd C. de Vries
    Mar 12, 2012 at 17:43
  • 2
    $\begingroup$ @SjoerdC.deVries Actually, it does prevent the function from being evaluated outside the region. Try for example something like f[x_, y_] := x^2 + y; Reap[NIntegrate[Boole[x + y < .5] f[x, y], {x, 0, 1}, {y, 0, 1}, EvaluationMonitor :> Sow[{x, y}]], _, ListPlot[#2] &], then the plot shows that the integrand is only evaluated in the region where Boole[]!=0. $\endgroup$
    – Heike
    Mar 12, 2012 at 18:01
  • $\begingroup$ @Heike OK, I didn't check that. I based my conclusion on the timing. In the given example here adding Boole really doesn't work to reduce calculation time. $\endgroup$
    – Sjoerd C. de Vries
    Mar 12, 2012 at 18:06
  • $\begingroup$ @SjoerdC.deVries I guess it depends how expensive the calculation of myF is whether it's worthwhile to use Boole or not. $\endgroup$
    – Heike
    Mar 12, 2012 at 18:11
  • 2
    $\begingroup$ @SjoerdC.deVries I have a hunch that when you add Boole, it's the symbolic preprocessing that takes time. Mathematica does a lot of smart things including checking the symmetries of the domain and the integrand to reduce the number of points where it needs to be evaluated. Check the link I edited in (I was surprised). Constraining the domain also influences the location of points where its evaluated (i.e. it doesn't just numerically filter out the points that are outside the domain---it seems to use symbolic methods to generate them inside the domain.) $\endgroup$
    – Szabolcs
    Mar 13, 2012 at 5:26
11
$\begingroup$

For total control over integration, sweep over the ellipse along the eigendirections. To get the integral right, it is essential to use unit-length eigenvectors. This code shows a discrete version of the sweep, to help you visualize it:

m = {{2, 1}, {1, 3}}; (* Matrix *)
c = {3, 5};           (* Center *)

ev = Sqrt[Eigenvalues[m]]; (* Semi-axes *)
evec = Eigenvectors[m];    (* Principal directions, unnormalized *)
evec = DiagonalMatrix[1/Norm[#] & /@ evec].evec; (* Unit eigenvectors *)
With[{k = 20}, 
 ListPlot[Flatten[
   Table[c + {s, t}.evec,
     {s, -ev[[1]], ev[[1]], ev[[1]]/k}, 
     {t, -ev[[2]] Sqrt[1 - s^2 / ev[[1]]^2], ev[[2]] Sqrt[1 - s^2 / ev[[1]]^2], 
          ev[[2]]/k }
   ], 1], 
  AspectRatio -> Sqrt[ev[[1]] / ev[[2]]], AxesOrigin -> {1.5, 3}, 
  PlotStyle -> Gray, 
  Epilog -> {PointSize[0.02], Point[c], Thick, Darker[Red], 
    Arrow[{c, c + {ev[[1]], 0}.evec}], Darker[Blue], 
    Arrow[{c, c + {0, ev[[2]]}.evec}]}]
 ]

This draws the first axis in red and the second in blue, each to the tip of the ellipse, on top of the dots created by the double-sweep across the eigendirections:

Ellipse

Integration is done in a similar fashion. Here's the area of the ellipse:

cnst = NIntegrate[1, {s, -ev[[1]], ev[[1]]}, 
         {t, -ev[[2]] Sqrt[1 - s^2/ev[[1]]^2], ev[[2]] Sqrt[1 - s^2/ev[[1]]^2]}]

It returns 7.02481, equal to $\sqrt{5}\pi$, as it should (the determinant of m is $5$).

As a further check we can recover the ellipse parameters from its first two moments. E.g.,

Clear[f];
f[x_, {p_, q_}] := x[[1]]^p x[[2]]^q; (* For the (p,q) moment *)
ParallelTable[
 NIntegrate[
   f[c + {s, t}.evec, pq], 
     {s, -ev[[1]], ev[[1]]},
     {t, -ev[[2]] Sqrt[1 - s^2 / ev[[1]]^2], ev[[2]] Sqrt[1 - s^2 / ev[[1]]^2] }
  ] / cnst, 
  {pq, {{1, 0}, {0, 1}}}]

This returns the center, {3., 5.}, as it should.

ParallelTable[
 4 NIntegrate[
    f[{s, t}.evec, pq], 
      {s, -ev[[1]], ev[[1]]},
      {t, -ev[[2]] Sqrt[1 - s^2 / ev[[1]]^2], ev[[2]] Sqrt[1 - s^2 / ev[[1]]^2] }
   ] / cnst, 
   {pq, {{2, 0}, {1, 1}, {0, 2}}}]

These three central second moments (scaled up by 4) ought to return the coefficients of the matrix m. Indeed, the result is {2., 1., 3.}, exactly as expected. The visualization and these checks should give you confidence that this technique is correct.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks: I fixed the expression and then checked it by pasting it from here back into Mathematica and executing it. BTW, replacing NIntegrate by Integrate works in all these examples (but takes a few seconds). $\endgroup$
    – whuber
    Mar 13, 2012 at 14:43
  • 1
    $\begingroup$ Why don't you use Orthogonalize to normalize the eigenvectors? As they're guaranteed to be orthogonal already, it will just normalize them. Also, there's Normalize which can be mapped across them to do the same thing. $\endgroup$
    – rcollyer
    Mar 13, 2012 at 14:56
  • 1
    $\begingroup$ Thanks for the suggestions--I am grateful to learn about Normalize and welcome any code improvements you might ever have the inclination to volunteer in the future. (It seems Mathematica has loads of stuff built in, but first you have to discover it.) But in this case we probably shouldn't care about elegance or efficiency, because computing the integration directions and limits is O(1) overhead for what promises to be a long painful integration; the primary concern therefore has to be correctness of code. (If it weren't for that, this answer would have been about two lines long. :-) $\endgroup$
    – whuber
    Mar 13, 2012 at 17:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.