Integration in 2D (area) using the Monte Carlo method
18 次查看(过去 30 天)
显示 更早的评论
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.
0 个评论
采纳的回答
Roger Stafford
2014-2-23
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.
更多回答(1 个)
Youssef Khmou
2014-2-23
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
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!