MATLAB: Integration in 2D (area) using the Monte Carlo method

circleellipsemonte carlo methodnumerical integration

I need to numerically compute the area of the region between an ellipse (say with major axis a and minor axis b) centered at origin and a concentric circle of radius R. Area of the region between the circle and the ellipse is clearly,

\pi R^2 - \pi ab

I realized that using the usual Gauss-Legendre quadrature for this does not give correct results. I divided the region into 4 sectors and tried to compute the area by using the GL quadrature and I end up with more than 2% error. The reason is the integration points are not correctly distributed and this happens because of the different radii of the ellipse and the circle (see picture). The result of this is a small area remaining unaccounted for and causing the error in the final result (one can see the white-space at 45 degrees near the ellipse boundary). For this figure, a 30 x 30 Gauss rule is used in the sector in first quadrant and the red dots are the Gauss integration points (the points that are seen clustering towards the end of the ellipse do not actually cross over into the ellipse). Is it possible to use the Monte Carlo method for computing the integral in this case? If yes, could someone point me to a working example in Matlab for computing a 2D integral. There are probably other elegant methods to correctly evaluate the area of the annular region between ellipse and the circle like conformal mapping but I think they are more challenging to implement. Thanks for the help.

Best Answer

Using a Monte Carlo method to solve an area like this would be a terrible method to use on that problem, either requiring huge numbers of random points or being fraught with relatively large errors.

Furthermore it seems a shame to even use numerical quadrature when the integrals involved have an explicit solution which will yield the answer you quoted, pi*(R^2-a*b), exactly.

But if you are intent on doing it numerically, then it should be set up correctly. If, as it appears in your image, you are doing it using polar coordinates, r and phi, the inner radial limit should have been

   a*b/sqrt(b^2*cos(phi)^2+a^2*sin(phi)^2)

I have no idea what you were using as an inner limit here but it is very wrong, and that is not the fault of the Gauss-Legendre integration method. Also it is important that you use the correct Jacobian factor, which with polar coordinates would be just r.

It is not essential that you use Gauss-Legendre integration here. Since the circle and ellipse are entities for which precise values can easily be calculated, there is no great disadvantage in using larger numbers of evenly-spaced points as in trapezoidal integration. Gauss-Legendre is most useful when there is a heavy computation penalty to pay for calculating an integrand and the number of such computations needs to be kept to a minimum.

Also note that it is silly to use numerical integration along the radial lines for the inner integral when you are merely computing the integral of the Jacobian, r. Its indefinite integral is r^2/2 so the definite integral would be the difference between r^2/2 at the outer limit, R, and the inner limit which I gave above.

In computing the resulting outer integral with respect to phi you might need to use numerical methods. However, if you had used cartesian coordinates originally instead of polar coordinates, the whole problem would have been easier and numerical methods would not be needed at all.

Related Solutions

MATLAB: Integration over a circle VS rectangle

Dont forget the basics:

dxdy = rdrdtheta

HTH

MATLAB: Numerically integrate function over 2D ellipse

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.

Best Answer

Related Solutions

MATLAB: Integration over a circle VS rectangle

MATLAB: Numerically integrate function over 2D ellipse

Related Question