MATLAB: Numerically integrate function over 2D ellipse

MATLABnumerical integrationsurface integral

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!

Best Answer

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.
Related Question