Numerically integrate function over 2D ellipse?

I am simulating a sensor that is assumed to sample concentration values C(x,y) in an arbitrary 2D elliptical area given by

 ( (x-x_s)*cos(alpha) - (y-y_s)*sin(alpha) )^2/a^2 + ( (x-x_s)*sin(alpha) + (y-y_s)*cos(alpha) )^2/b^2

where (x_s, y_s) is the center of the ellipse,

alpha is the rotation of the ellipse, and

a and b are the axes.

I need to estimate the surface integral of C(x,y) over this ellipse.

An old post on Loren on the Art of Matlab gives some hints for setting up numerical integration over an ellipsoid region, but I am unsure how to proceed. I have also looked at the integral2 function, which was introduced after Loren's blog post, but am unsure how to define the x,y bounds to account for the ellipse.

For the sake of testing, I define C(x,y) as a simple function:

conc = @(x,y) exp(-sqrt(x.^2 + y.^2) /2)

The integral2 function behaves as expected when I give it rectangular bounds:

integral2(conc,-1,1,-1,1)

ans =

    2.7565

Can anyone provide advice on how to modify the domain of interest in integral2? Or provide an alternative method? Thanks!

Answer 1

[Ameer Hamza]

If you look at the first example of the documentation page integral2(). You will get the idea, how to integrate over an ellipse. But the problem here is that you need an explicit equation of the ellipse in the form y=f(x). But in your case, the ellipse is rotated, therefore its equation cannot be expressed in term of y=f(x). You will need to find a way in which you can convert your integration function into a form where you will get an explicit equation of the ellipse. In the following answer, I am using the following values of the parameters

x_s = 1;
y_s = 2;
alpha = pi/6;
a = 2;
b = 1;

1) We need to remove the rotation of the ellipse. To remove rotation, we first need to translate the center of the ellipse to the origin. So we need to apply following two transformation

1) Translation: 
   x -> x+x_s  and y -> y+y_s  
This transformation will move the ellipse to the origin.

2) Rotation: Rotate axis by alpha degree counter-clockwise. 
This will align the ellipse with x-y axis.

Your original ellipse is like this

After the first transformation (i.e. translation) it becomes like this

And after second transformation, it becomes

after these transformations, the equation of the ellipse becomes

x^2/a^2 + y^2/b^2 = 1

now this is not a rotated ellipse anymore and we can write an explicit equation of the ellipse as following

y = b*sqrt(1-x^2/a^2)

2) Since we applied the transformation to the ellipse, to conserve the value of integral, we must apply these two transformations to concentration values conc() also. Therefore the original surface looks like this

applying the first transformation (i.e. translation) the equation becomes

conc = @(x,y) exp(-sqrt((x+x_s).^2 + (y+y_s).^2) /2);

and the surface becomes

this is shifted as compared to the first surface. Now applying the rotation (alpha angle), the equation becomes

conc = @(x,y) exp(-sqrt((x*cos(alpha)-y*sin(alpha)+x_s).^2 + (x*sin(alpha)+y*cos(alpha)+y_s).^2) /2);

and the surface becomes

3) Now do numerical integration with the modified conc function. First define the ellipse

ellipse = @(x) b*sqrt(1-x.^2/a.^2);

but this is just the top half of ellipse, also define the bottom half of ellipse

ellipse_ = @(x) -b*sqrt(1-x.^2/a.^2);

Now use integral2 as follow

result = integral2(conc, -a, a, 0, ellipse) + integral2(conc, -a, a, ellipse_, 0)
result =
  2.1063

Here conc is the final expression we obtained, after translation and rotation.

[James Hengenius]

Thank you! I was so focused on the translation/rotation of the ellipse (sensor) in space, it didn't occur to me to translate the concentration function.

Your step-by-step explanation was clear and very helpful - thank you!

[Ameer Hamza]

You are welcome.

[John D'Errico]

+1. A reasonable answer.

Personally, I'd just view this as a conversion to polar coordinates. An ellipse is easy enough to write in a polar form, so r varies according to theta, thus polar angle. Now write the double integral over the ellipse in polar form. Don't forget the dx*dy turns into r*dr*dtheta. Integrate over polar angle, varying from 0 to 2*pi, and r from 0 to the perimeter of the ellipse.

A simple ellipse is of the form

x^2/a^2 + y^2/b^2 = 1

If we have x=r*cos(theta), and y=r*sin(theta), then we can simply solve for r as a function of theta. So

r^2 = 1/(cos(theta)^2/a^2 + sin(theta)^2/b^2)

A rotated, translated ellipse is just as easy to write. So the integral on r takes on limits of 0 to the positive sqrt of the right hand side.

Finally, note that if the ellipse is translated, then you simply define the center of the polar coordinate form as the center of the ellipse, translating the entire problem.

[James Hengenius]

Thank you! I was so focused on the translation/rotation of the ellipse (sensor) in space, it didn't occur to me to translate the concentration function.

Your step-by-step explanation was clear and very helpful - thank you!

John - thank you for noting the integral in polar coordinates. However, the spatial integral of my concentration function will not always have an analytic solution (e.g., interpolated from data, I should have mentioned this), which is why I was looking for a numerical solution such as integral2.

Answer 2

[John D'Errico]

James. You clearly have not thought about what I said. I never said you need to use an analytical integration. All I showed was how to simply transform the integration into polar coordinates, with the origin translated to the center of the ellipse. You will still use integral2!!

The point I made is that when you integrate in polar coordinates, then you will need an extra factor of r in the kernel.

Do I really need to explain this in great detail? Argh. Since I apparently need to show you how to do this in detail, I'll write it as an answer.

For example, integrate the function x*y, over an elliptical domain, centered at (x,y) = [1 2]. The equation of the ellipse will be a simple one:

(x-1)^2 + (y-2)^2/4 = 1

Thus the ellipse has a major axis that is twice the length of the minor axis.

xycenter = [1 2];
fxy = @(x,y) x.*y; % You supply this

fthetar = @(theta,r) fxy(xycenter(1) + r.*cos(theta),xycenter(2) + r.*sin(theta));

So nothing of any significance required there. All you need to supply is the basic function fxy, and the center of the ellipse.

I already gave you the formula for the ellipse boundary in polar form, at least for an unrotated ellipse, of the same form, but since I've already dealt with the translation, we can ignore that part. If the ellipse was rotated, we would have an xy term in the ellipse equation.

a = 1;
b = 2;
rbound = @(theta) sqrt(1./(cos(theta).^2/a^2 + sin(theta).^2/b^2));

Does it work? TRY IT!

rbound(0)
ans =
     1

rbound(pi/2)
ans =
     2

So at 0 degrees, we get a boundary in r of 1. The ellipse is elongated in the y direction, so at 90 degrees, it goes out to 2.

Now, integrate. USE INTEGRAL2. See that I carefully multiply by r in there.

I = integral2(@(theta,r) fthetar(theta,r).*r,0,2*pi,0,rbound)
I =
        12.5663706004049

Is that correct? Perhaps simplest is to use an alternative method. Here I'll show an even simpler solution, a tool called tessquad, that does a gaussian integration over a triangulated domain. Apparently it is not posted on the file exchange.

I'll create a set of points around the perimeter of the ellipse, then build a triangulated domain from that polygon.

nt = 100;
t = linspace(0,2*pi,nt)';
xyellipse = [xycenter(1) + rbound(t).*cos(t),xycenter(2) + rbound(t).*sin(t)];

% ALAYS PLOT EVERYTHING
plot(xyellipse(:,1),xyellipse(:,2),'-')
axis equal
grid on

Now, build a triangulation over that domain

vert = [xyellipse;xycenter];
tri = [(1:nt)',[(2:nt)';1],repmat(nt+1,nt,1)];

Plot. AGAIN.

trisurf(tri,vert(:,1),vert(:,2),fxy(vert(:,1),vert(:,2)))

That is a crude plot of the surface. Note that the actual integration will be exact, because it will employ a Gaussian integration of sufficient order, and this is a simple polynomial function on that domain. Since the function is such low order, the gaussian integration is a very basic one, but tessquad worries about that.

It = tessquad(@(xy) fxy(xy(:,1),xy(:,2)),1,vert,tri)
It =
         12.550841380856

Or, given the triangulation, I can do it using an adaptive 2-dimensional quadrature tool I've recently written to work on a triangulation.

options.tol = 1e-12
options = 
  struct with fields:
    tol: 1e-12

I = simplxint2(fxy,tri,vert,options)
I =
         12.550841380856

These two solutions have generated the exact solution to an integral over a triangulated domain that approximates the ellipse, because the boundary is represented as a polygon, not the true ellipse. Both of the above solutions are only approximations, because they miss a tiny bit of the volume. Had I made nt larger, they would be now much closer. I deliberately set nt at only 100 before so the surface plot would be best seen. Otherwise the edges in the triangulation would have dominated the surface.

nt = 1000;

(having redone the triangulation to reflect the larger nt)

It = tessquad(@(xy) fxy(xy(:,1),xy(:,2)),1,vert,tri)
It =
        12.5662178639005

Now, exactly where in here do you think I've done any analytical integration? All you needed to supply was the function fxy, the center of the ellipse, and the lengths of the ellipse axes. In this case, I used a basic, unrotated ellipse. But that is trivially changed.

I should repost tessquad and simplxint2 on the FEX I guess, but they were used here only to verify the result from integral2.