pythonmathsimulationcomputational-geometrymontecarlo

Calculate the area of ​intersection of curves using the Monte Carlo Method in Python


I'm making a list of exercises in which I have the following exercise:

"Use the Monte Carlo method to calculate the area of ​​intersection between the curves given below. For random number generation use the Linear Congruence Method"

curves functions

I'm having trouble determining the bounds of random numbers. I already added the functions in GeoGebra


Solution

  • Define box which contains whole area: I would say (-0.5,-0.5) to (1.5,2.0) should do it.

    Sample points in the box:

    x = -0.5 + 2*np.random.random()
    y = -0.5 + 2.5*np.random.random()
    

    Check against four curves if point is inside area

    def inside(x, y):
        if ... :
             return True
        return False
    

    Count ratio of points inside vs total points and multiply by box area

    N = 100000
    is_inside = 0
    
    for k in range(0, N):
        x = -0.5 + 2*np.random.random()
        y = -0.5 + 2.5*np.random.random()
    
        if inside(x,y):
            is_inside += 1
    
    area = is_inside*(2.5+0.5)*(1.5+0.5)/N
    

    UPDATE

    Anyway, this is good solution, I think. Python 3.10 x64 Windows 10

    import numpy as np
    import matplotlib.pyplot as plt
    
    rng = np.random.default_rng(135797531)
    
    def inside(x: np.float64, y: np.float64) -> bool:
    
        if x <= 0.0:
            return False
    
        y1 = 4.0*x
        if y > y1:
            return False
    
        y2 = x*np.log(x)
        if y < y2:
            return False
    
        y3 = np.cosh(x)
        if y > y3:
            return False
    
        y4 = np.cos(x) / x
        if y > y4:
            return False
    
        return True
    
    x_min: np.float64 =  0.0
    y_min: np.float64 = -0.5
    
    x_max: np.float64 = 1.5
    y_max: np.float64 = 2.0
    
    N:int = 10000
    
    is_inside:int = 0
    
    x_plt = []
    y_plt = []
    
    for k in range(0, N):
        
        x = x_min + (x_max - x_min)*rng.random()
        y = y_min + (y_max - y_min)*rng.random()
    
        if inside(x,y):
            is_inside += 1
            x_plt.append(x)
            y_plt.append(y)
        
    area = (x_max - x_min)*(y_max - y_min)*is_inside/N
    print(area)
    
    plt.scatter(x_plt, y_plt, s=1)
    plt.show()
    

    that will produce output equal to 1.18275 and plot below

    enter image description here