2
$\begingroup$

Suppose I want to use the Monte Carlo integration method to compute the following integral

$\int_D (e^{x^{2}} + e^{y^{2}}) \; dx \; dy$ where $D$ is some regular hexagon. I have managed to write the code and everything and the results I'm getting seem valid.

The problem is now that I do not know how can I show that the algorithm does what it's supposed to do. Is there another way to compute that integral? Is there any way to validate my results?

EDIT: I'm talking about the "sample-mean" method

  • 0
    Perhaps see `http://math.stackexchange.com/questions/1635250/` for more on Monte Carlo integration.2016-02-01

3 Answers 3

0

Because the integrand is positive you can easily find lower and upper bounds by integrating over nicer domains that are subsets or supersets of the hexagon domain.

  • 0
    While this is true, such an approach require fine mesh triangulation of the hexagon to achieve modest precision, as $\exp(x^2) + \exp(y^2)$ grow rather fast.2012-01-25
1

The simplest way to validate your Monte Carlo results would be to use numerical quadrature of modestly accurate order on the region D. Since D is a regular hexagon (but otherwise unspecified), it can be split into six equilateral triangles.

Quadrature rules in one-dimension are well introduced in undergraduate numerical methods courses, but if you want, take a formula from this article on quadrature rules for triangles.

Typically a series of such approximations is done, either with increasingly small subdivisions of area or with increasingly higher-order accuracy formulas, so that the convergence of the approximations can be seen.

1

Let $f(x,y) = \mathrm{e}^{x^2} + \mathrm{e}^{y^2}$, and suppose we would like to integrate this function over a valid triangle with vertices at $p_1 = (x_1,y_1)$, $p_2 (x_2,y_2)$ and $p_3 = (x_3,y_3)$.

This can be computed in closed form, because $\int_D f(x,y) \mathrm{d} x \mathrm{d} y = \int_D \mathrm{e}^{x^2} \mathrm{d} y \mathrm{d} x + \int_D \mathrm{e}^{y^2} \mathrm{d} x \mathrm{d} y $ Now integration over $y$ for each fixed $x$ in the first integral, and over $x$ for each fixed $y$ in the second can be easily done, and will yield a linear function of $x$ and $y$ respectively.

Then you are down to $ \int_a^b (c x+d) \mathrm{e}^{x^2} \mathrm{d} x = \frac{c}{2} \left( \mathrm{e}^{b^2}- \mathrm{e}^{a^2} \right) + \frac{d \sqrt{\pi}}{2} \left( \operatorname{erfi}(b) - \operatorname{erfi}(a) \right) $

Coming back to the setting stated at the beginning of the post, here is a little code in Mathematica that will do the computation exactly:

In[6]:= pred[{x_, y_}, {x1_, y1_}, {x2_, y2_}, {x3_,     y3_}] := ((x - x1) (y2 - y1) - (y - y1) (x2 - x1)) ((x3 - x1) (y2 -          y1) - (y3 - y1) (x2 - x1)) > 0  In[7]:= InTriangle[px_, p1_, p2_, p3_] :=   pred[px, p1, p2, p3] && pred[px, p2, p3, p1] && pred[px, p3, p1, p2]  In[10]:= int[p1 : {x1_, y1_}, p2 : {x2_, y2_}, p3 : {x3_, y3_}] :=   Block[{x, y},   Integrate[(Exp[x^2] + Exp[y^2]) Boole[      InTriangle[{x, y}, p1, p2, p3]], {x, Min[x1, x2, x3],      Max[x1, x2, x3]}, {y, Min[y1, y2, y3], Max[y1, y2, y3]}]]  In[13]:= int[{0, 0}, {1, 0}, {1/2, Sqrt[3]/2}] // FunctionExpand  Out[13]= 1/6 (-Sqrt[3] + 6 Sqrt[3] E^(1/4) - 2 Sqrt[3] E^(3/4) -     3 Sqrt[3] E - 3 Sqrt[3 Pi ] Erfi[1/2] +     3 Sqrt[3 Pi] Erfi[1] + 3 Sqrt[ Pi] Erfi[Sqrt[3]/2])