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

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.

Answer 1

[Roger Stafford]

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.

[Ganesh]

Thanks for your comments.

_ It is not essential that you use Gauss-Legendre integration here. _

I agree, but, I need to use Gauss-Legendre scheme as I am using the classical Finite Element shape functions that are defined on [-1, 1] intervals. Although, this is not a necessity, but using different shape functions would mean making significant changes in my main code.

_ I have no idea what you were using as an inner limit here but it is very wrong _

Thanks for pointing this out. I made the changes (see the attached code) and the area computed now matches reasonably well with the analytical value. The error in computation now is 2.23E-07%. But, the integration points still do not appear to be correctly placed. This means the limits on the radius is still wrong! I am clearly doing some really basic mistake here. Appreciate your time and help.

[Ganesh]

I had made a mistake in the code. The integration points now appear OK and the percentage error in the area is 5.9E-014. All good. Thanks a lot for pointing out my mistake earlier. The new code is attached in case anyone is interested.

Answer 2

[Youssef Khmou]

Ganesh,

In this second answer to the problem, i wrote a Monte Carlo test for this specific geometry , try to compute the surface with and answer me back :

% Monte Carlo method for computing the area between
% an ellipse (a,b) and quart of cirle (r>(a,b)),  initiation 
clear;
r=5;
a=2;
b=1;
% x/a^2+y/b^2=1
N=10000;
A=0;
for i=1:N
  p=r*abs(randn(1,2));
  x=p(1);y=p(2);
  q1=((x/a).^2)+((y/b).^2);
  q2=sqrt((x.^2)+(y.^2));
  if (q1>=1.00 && q2<=r)
      A=A+1;
      plot(x,y,'*');
      hold on;
  end
end

[Ganesh]

Thanks for your answer. My apologies but as I mentioned in my reply to Roger, I may not be able to move away from the Gauss-Legendre quadrature because of my FE shape functions which are defined on a canonical square [-1, 1] X [-1, 1]. I should confess that I do not understand how Monte Carlo method works but should I use it, I think I will still need to map the random points on to a [-1, 1] X [-1, 1] space. I will then need to compute the Jacobian at each random point and then finally compute the area.