Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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

share|improve this question
This is just a guess, but maybe you have a typo and you really want to integrate over $e^{-(x^2 + y^2)}$ instead? This would be a 2d gaussian which is probably more common than the function you wrote. –  opt Jan 25 '12 at 18:54

3 Answers 3

up vote 1 down vote accepted

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])
share|improve this answer

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.

share|improve this answer
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. –  Sasha Jan 25 '12 at 18:40

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.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.