...

Batman Integration

Background

The Monte Carlo algorithm is used all over the world for a number of things. In physics, it is often used to find integrals of poorly-behaved functions, and it works incredibly well for that! Let’s take a minute to describe how the Monte Carlo algorithm works and why it’s cool.

Firstly, let’s define what an integral is: it is a fancy way to find an area. For now, that’s all we need to know. Sure, there are plenty of people who make it way, way more complicated than that, but the simpler the definition, the simpler the algorithm. I like simple algorithms.

So, if an integral is just an area, what do we need to do to find that area? Well, use geometry! If we have a box, the area of the box is obviously its . If we have a triangle, the area is simply . Generally, the more complicated the shape, the harder it is to find the shape’s area from geometry.

That said, we don’t really need to find the shape of the object we are integrating to find it’s area. Instead, we can take a simple object, find the area of that and then embed a more complicated object inside of it. If we can find the ratio of how filled the simple object is by the more complicated object inside of it, then we can find the area of the complicated object like so:

Of course, this will only work if the simple object is larger than the more complicated one, such that the complicated object falls perfectly inside of the larger object’s boundaries.

Let’s take the example of a circle and square:

Circle inside of Square

In this case, the area of the circle is and the area of the square is . This means that

From here, it’s obvious that

Or the area of the circle.

But now we have a question: How do we find that ratio for arbitrary shapes? Well, Monte Carlo’s solution to this is Random Sampling. If we take a bunch of points in the simple region with the more complicated shape inside, then

And that’s pretty cool. It means that if we take a whole bunch of points, we should be able to find the ratio and integrate any arbitrarily shaped / size object! For the example above:

gif of simple Monte Carlo algorithm

And we see that the ratio is pretty darn close to .

The most difficult part of the Monte Carlo algorithm is finding the appropriate evaluation function. For example: How did we know whether a random point on the plane was within the circle above? Well, we simply had to plug it into a function that checked to see whether a line drawn from the origin to the x-y location of the random point was greater or less than a provided radius for the circle. The code basically looks like:

// Function to check whether point is in circle or not
bool in_circle(pos loc, pos origin, double radius){

    double x = loc.x - origin.x;
    double y = loc.y - origin.y;

    if (x*x + y*y < radius * radius){
        return true;
    }
    else{
        return false;
    }
}

That said, the true power of the Monte Carlo algorithm comes from the fact that we may now integrate somewhat arbitrary functions, for a good example of this, see the next section!

Integrating the Batman Curve

Now here’s the thing: Monte Carlo doesn’t need to use boxes and it doesn’t need to integrate obvious shapes. Let’s talk about what you need to do to integrate non-obvious functions with one of our favorite superheroes: Batman!

To be clear: the Batman Curve was not created by me. I found it on reddit and found all the information I needed on this math stack exchange post. In addition, WolframMathWorld had a decent page with the integral on it. Because the Batman Curve is difficult to integrate (though not impossible by any means!) and it has a clear analytically-defined integral to compare our Monte Carlo integration algorithm against, it was the perfect choice to test some things out.

First things first, for this integral, we were not using a standard box as described above. Instead we were using an oval. Honestly, the objects don’t matter, but I thought it was worthwhile to point out that the area of the oval is simple , so with the shortest and longest distances to the oval from the origin as the radii.

The batman curve was a little tricky, we had to split it into the following equations (image created by Jack in the previously mentioned math stack exchange post):

Image of equations, found here: http://math.stackexchange.com/questions/54506/is-this-batman-equation-for-real

Which, when plotted provide:

Batman Curve

Then we had to write a simple function to determine if a point inside the oval was within the Batman Curve or not. To do this, we divided the curve into the following sections:

Batman Curve with blocks

And here is a visualization of the output:

gif of Batman Monte Carlo

Here is the function (Nothing special, really):


// Function to check if we are in batman function for monte_carlo
bool is_batman(pos dot, pos ori, double scale){

    // Determine x and y locations
    double pos_x = (dot.x - ori.x) / scale;
    double pos_y = (dot.y - ori.y) / scale;
    double temp_x, temp_y;

        // Drawing left upper wing
        if (pos_x <= -3){
            temp_x = (-7 * sqrt(1-((pos_y * pos_y)/9.0)));
            if (pos_x >= temp_x){
                return true;
            }
            else{
                return false;
            }
        }

        // Drawing left shoulder
        if (pos_x > -3 && pos_x <= -1){
            temp_x = -pos_x;
            temp_y = -(2.710523708715754 + (1.5 - 0.5 * temp_x)) 
                     +1.355261854357877
                     *sqrt(4.0-(abs(temp_x)-1)*(abs(temp_x)-1));
            if (pos_y > temp_y){
                return true;
            }
            else{
                return false;
            }
        }

        // Drawing exterior left ear
        if (pos_x > -1 && pos_x <= -0.75){
            temp_y = 9.0 + 8 * pos_x;
            if (pos_y > -temp_y){
                return true;
            }
            else{
                return false;
            }
        }

        // Drawing interior left ear
        if (pos_x > -0.75 && pos_x <= -0.5){
            temp_y = -3 * pos_x + 0.75;
            if (pos_y > -temp_y){
                return true;
            }
            else{
                return false;
            }

        }

        // Drawing top of head
        if (pos_x > -.5 && pos_x <= 0.5){
            temp_y = 2.25;
            if (pos_y > -temp_y){
                return true;
            }
            else{
                return false;
            }

        }


        // Drawing interior right ear
        if (pos_x > 0.5 && pos_x <= 0.75){
            temp_y = 3 * pos_x + 0.75;
            if (pos_y > -temp_y){
                return true;
            }
            else{
                return false;
            }

        }

        // Drawing exterior right ear
        if (pos_x > 0.75 && pos_x <= 1){
            temp_y = 9.0 - 8 * pos_x;
            if (pos_y > -temp_y){
                return true;
            }
            else{
                return false;
            }
        }


        // Drawing right shoulder
        if (pos_x <= 3 && pos_x > 1){
            temp_y = -(2.710523708715754 + (1.5 - 0.5 * pos_x)) 
                     +1.355261854357877
                     *sqrt(4.0-(abs(pos_x)-1)*(abs(pos_x)-1));
            if (pos_y > temp_y){
                return true;
            }
            else{
                return false;
            }
        }

        
        // Drawing right upper wing
        if (pos_x > 3){
            temp_x = (7 * sqrt(1-((pos_y * pos_y)/9.0)));
            if (pos_x <= temp_x){
                return true;
            }
            else{
                return false;
            }

        }
    }
    if (pos_y >= 0){
         // Drawing bottom left wing
         if (pos_x <= -4.0){
             temp_x = (-7 * sqrt(1-((pos_y * pos_y)/9.0)));
             if (pos_x >= temp_x){
                 return true;
             }
             else{
                 return false;
             }
         }
         // Drawing bottom wing
         if (pos_x > -4.0 && pos_x <= 4){
             //return false;
             temp_y = (abs(pos_x/2) 
                       - 0.09137221374655434 * pos_x * pos_x - 3.0)
                      + sqrt(1 - (abs(abs(pos_x)-2)-1) 
                             * (abs(abs(pos_x)-2)-1));
             temp_y *= -1;
             if (pos_y < temp_y){
                 return true;
             }
             else{
                 return false;
             }
         }
         // Drawing bottom right wing
         if (pos_x >= 4.0){
             //return true;
             temp_x = (7 * sqrt(1-((pos_y * pos_y)/9.0)));
             if (pos_x <= temp_x){
                 return true;
             }
             else{
                 return false;
             }
         }
    }

    return false;
    
}

The full code is available here. Note that it is a little messy and long. That’s because the we used Cairo to do the visualization, which added a bunch of extra lines to the code.

Overall, this was a pretty fun side-project!

Thanks for reading!