Numerically integrate function over 2D ellipse?

16 次查看（过去 30 天）
James Hengenius
James Hengenius 2018-5-21
编辑： John D'Errico 2018-5-29
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!

请先登录，再回答此问题。

采纳的回答

Ameer Hamza
Ameer Hamza 2018-5-21
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.
  4 个评论
显示 2更早的评论
John D'Errico
John D'Errico 2018-5-25
编辑：John D'Errico 2018-5-25
+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
James Hengenius 2018-5-29
编辑：James Hengenius 2018-5-29
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.

请先登录，再进行评论。

更多回答（1 个）

John D'Errico
John D'Errico 2018-5-29
编辑：John D'Errico 2018-5-29
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.

类别

在 Help Center 和 File Exchange 中查找有关 Transmitters and Receivers 的更多信息

标签

产品


版本

R2018a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by