Monte Carlo Simulations: How Big Is Your Heart?

Our final problem has no known exact solution. We want to find the area of the shape formed by,

(x^{2} + y^{2} - r^{2})^{3} - a\,x^{2}y^{3} \leq 0

The inequality has two parameters $r$ and $a$ . They are quantities of length, so they take on a nonnegative value. Let's try out some values and see what these parameters do. I have colored the resulting shapes in red. (Note, the scales are different.)

When $r = a$ (move along the diagonal, starting from bottom-left), we get the shape of a nice heart. When we increase $r$ while holding $a$ fixed, the heart morphs into a circle. On the other hand, when we increase $a$ while holding $r$ fixed, the heart turns into two petals. Well, I see bunny ears.

Clearly, the shape (i.e. area) of our heart depends on the radius $r$ and the ear length $a$ . Can you guess the formula for the area $A = A(r, a)$ ?

1. Problem description

Let $r,\,a \geq 0$ . Write down the expression for the area of the heart, the region that satisfies the inequality,
$(x^2 + y^2 - r^2)^3 - a\,x^2y^3 \leq 0$
The parameters $r$ and $a$ represent the radius and ear length of the heart.

I find this problem very interesting because it looks elegant and simple, yet we cannot (at least, I can't) find the area analytically.

Try writing the equation (the boundary of the heart) in polar coordinates $(R,\,\theta)$ :

(R^2 - r^2)^3 - a\,R^5\cos^{2}(\theta)\sin^{3}(\theta) = 0

If we can write $R$ as a function of $\theta$ , then we can apply the area formula

\displaystyle A(r,\,a) = \int_{\theta_{1}}^{\theta_{2}}\,\frac{1}{2}\,R(\theta)^{2}\,\,d\theta

But solving the sextic equation for $R$ , or even $R^{2}$ , looks impossible.

2. Shape of my heart (cf. Sting)

Although we cannot find the area analytically, we can still say a few things about what the expression must look like. We do so using dimensional analysis and boundary value analysis.

a. Dimensional analysis

We know that an area has the dimension of $\text{(length)}^{2}$ . In the inequality, the only parameters with dimension of $\text{(length)}$ are the radius $r$ and the ear length $a$ .

Hence, the area of our heart must consist of powers of $r$ and $a$ , and the powers must add to 2:

\left\{\begin{array}{l} A(r,\,a) = C_{1}\,r^{m_{1}}a^{n_{1}} + C_{2}\,r^{m_{2}}a^{n_{2}} + C_{3}\,r^{m_{3}}a^{n_{3}} + \cdots \\ C_{i}\text{'s are dimensionless constants} \\ m_{i} + n_{i} = 2,\,\,\text{for all}\,\,i \end{array}\right.

I have written the area as a linear combination of power terms $r^{m_{i}}a^{n_{i}}$ , because there is more than one answer to finding powers $m$ and $n$ that add to 2. We could have "nice" terms like $r^{2}a^{0}$ , $r^{1}a^{1}$ , and $r^{0}a^{2}$ , terms with a negative power like $r^{-1}a^{3}$ , and terms with a decimal power like $r^{1.4}a^{0.6}$ . All of these have dimension of $\text{(length)}^{2}$ , so we can linearly add them to obtain the most general expression for the area.

(As an aside, decimal powers involve real numbers, of which there are uncountably many, so it is incorrect to write integer indices as I did out of convenience.)

The expression contains infinitely many terms, so it seems that we are still doomed. We now argue that only 3 of them are, in fact, possible.

First, terms with a negative power are not allowed. Soon we will show that the area of our heart makes sense even when the radius and ear length are 0. If the expression included terms with a negative power, the area would approach to infinity as the radius or ear length tends to 0. We did not see such outrageous behavior in our diagram at the beginning of this post.

Second, terms with a decimal power are not allowed. Look back at our diagram. Doesn't the shape of the heart transition smoothly as we vary the values of $r$ and $a$ ? We posit that, not only is the area bounded and continuous in $r$ and $a$ , but also its first partial derivatives. Terms with a fractional power result in an unbounded derivative as the radius or ear length tends to 0. We also can't have a decimal power greater than 2, since we already eliminated terms with a negative power.

In conclusion, the area consists of only three terms:

\boxed{A(r,\,a) = C_{1}\,r^{2} + C_{2}\,ra + C_{3}\,a^{2}}

There are three dimensionless coefficients $C_{1}$ , $C_{2}$ , and $C_{3}$ for us to find.

b. Boundary value analysis

Recall the equation of our heart and the expression for the area inside:

\left\{\begin{array}{l} (x^{2} + y^{2} - r^{2})^{3} - a\,x^{2}y^{3} = 0 \\ A(r,\,a) = C_{1}\,r^{2} + C_{2}\,ra + C_{3}\,a^{2} \end{array}\right.

A boundary value analysis considers the extremes of the problem. We can imagine two extreme cases.

What happens when $a = 0$ ?
What happens when $r = 0$ ?

When $a = 0$ , the heart equation simplifies to the equation of a circle,

x^{2} + y^{2} = r^{2}

and we know the area of a circle! Solve $A(r,\,0) = \pi r^{2}$ to get $C_{1} =\pi$ .

Next, consider $r = 0$ and write the heart equation in polar coordinates:

R^{5} \cdot \Bigl[R - a\,\cos^{2}(\theta)\sin^{3}(\theta)\Bigr] = 0

When $R = 0$ , the equation holds trivially and we do not know how $R$ is related to $\theta$ any more than we did before.

So consider $R > 0$ . We can divide both sides of the equation by $R^{5}$ and say that,

R = a\,\cos^{2}(\theta)\sin^{3}(\theta),\,\,\text{for}\,\,\theta \in [0,\,\pi]

I have restricted the angle $\theta$ so that the polar curve draws the two ears only once. We can now use the area formula:

\displaystyle \begin{aligned} A(0,\,a) &= \int_{0}^{\pi}\,\frac{1}{2}\,\,R(\theta)^{2}\,\,d\theta \\ &= \frac{1}{2}\,a^{2}\,\int_{0}^{\pi}\,\,\cos^{4}(\theta)\sin^{6}(\theta)\,\,d\theta \\ &= \frac{3\pi}{512}\,a^{2} \end{aligned}

Now we know $C_{3}$ .

In summary, by analyzing the two boundary cases, we were able to determine two unknown coefficients:

\boxed{\displaystyle A(r,\,a) = \pi r^{2} + C_{2}\,ra + \frac{3\pi}{512}\,a^{2}}

Only one coefficient remains unknown, $C_{2}$ . We use Monte Carlo simulation to estimate its value.

3. Monte Carlo simulation

We know the coefficient $C_{2}$ once we know the area of our heart:

\displaystyle C_{2} = \frac{A(r,\,a) - \pi r^{2} - \frac{3\pi}{512}\,a^{2}}{ra}

We estimate the area by randomly generating points inside a box $[-L,\,L] \times [-L,\,L]$ . The probability that a random point lies inside the heart (i.e. it satisfies the inequality) depends on how large the heart is relative to the box:

\displaystyle p = \frac{\text{area of the heart}}{\text{area of the box}} = \frac{A(r,\,a)}{(2L)^{2}}

If we consider a point inside the heart to be a "win," we can solve the equation for the area of the heart:

\displaystyle A(r,\,a) \approx (2L)^{2} \cdot \frac{\text{number of wins}}{N}

Note that $L$ must be large enough to fully contain the heart inside.

By trial and error, I found that we should take $L$ to be the maximum of $1.5r$ and $0.25a$ . Our 3-term model says that the shape of the heart changes linearly with respect to either the radius or ear length. The fact that $L$ also changes linearly is further proof that our model is valid.

The code below shows how to calculate the area of our heart. Once again, we use an affine transformation to generate points inside the box, and vectorization to compute the area fast.

function A = calculate_area(r, a, N)
  % Set the size of the box
  L = max(1.5*r, 0.25*a);
 
  % Generate N points in the box
  x = (2*L)*rand(1, N) - L;
  y = (2*L)*rand(1, N) - L;
 
  % Check if a point is inside the heart
  criterion = ((x.^2 + y.^2 - r^2).^3 - a * x.^2 .* y.^3 <= 0);
 
  % Count how many points are inside the heart
  numWins = sum(criterion);
 
  % Return the area of the heart
  A = (2*L)^2 * (numWins / N);
end

I used $r = a = 1$ and $N = 10^{7}$ , and found $C_{2} \approx 0.501$ . In summary, the area of our heart is given by,

\boxed{\displaystyle A(r,\,a) \approx \pi r^{2} + 0.501\,ra + \frac{3\pi}{512}\,a^{2}}

4. How good is our model?

Let's see how well our model predicts the area. We perform Monte Carlo simulation for various values of $r$ and $a$ , and take the resulting areas to be the truth. They are as close as we can get to the exact values. We compare the areas predicted by our model to those produced by Monte Carlo simulation.

With $N = 10^{4}$ simulations, we get the following:

We see that the area predicted by our model coincides with that from Monte Carlo simulation well, both in terms of the graph and the absolute error.

The absolute error, which lingers in the hundreds (sometimes over 1,000), seems large. However, relative to the area, it is not. When the radius is 100, Monte Carlo says that the area is around 35,000. Compared to this, an absolute error of 1,000 is quite small (2.86% relative error).

We also notice that the results from Monte Carlo simulation are less consistent when the radius is large. We can counteract the size effect by running more simulations for large radii and (to be safe) large ear lengths. That is, we allow $N$ to be a function of the parameters $r$ and $a$ :

\left\{\begin{array}{l} \alpha = \ln(1 + ra) \\ N(r,\,a) = \text{floor}\bigl(10^{4} \cdot (1 + \alpha)\bigr) \end{array}\right.

The number $\alpha$ tells us the factor by which we increase the base number of simulations. I have designed $\alpha$ so that we do not need additional simulations when the radius or ear length is 0, and gradually increase the number of additional simulations as we increase the radius or ear length.

With this simple change, we get better looking results:

Area produced by Monte Carlo simulation (variable N)

Absolute error between the two areas (variable N)

5. Conclusion

Over the last four posts, we studied what Monte Carlo simulation is and how we can use it to solve various problems. It's an extremely versatile tool, in that, we can adapt it to many different fields (we considered probability theory and integration, in particular) and can even use it to solve problems for which there is no exact solution.

Along the way, we also considered several Matlab's built-in functions and explored a few programming concepts, including control statements, user-defined functions, and polymorphism.

Our journey in understanding math has only just started. Until next time!

Notes

You can find the code in its entirety and in other languages here: Download from GitHub