[Scientific / technical calculation by Python] Monte Carlo integration, numerical calculation, numpy

Introduction

By Monte Carlo integration method Calculate $ I = \ int_0 ^ 1 \ frac {4} {1 + x ^ 2} dx = \ pi $.

Contents

(1) ** Calculation by simple Monte Carlo method **. The area is evaluated by how many generated random numbers enter the integration area (how many balls are thrown into the mat).

(2) ** Monte Carlo integration calculation by the mean method **. A method in which the domain of a function is sampled uniformly and the average value of y values at the generated x points is used. Generally, if the mean value of y is $ y_ {av} $,

Integral $ I = \ int_a ^ b f (x) dx \ sim \ frac {b-a} {N} y_ {av} $

Will be. Therefore, it is the miso of this method to determine $ y_ {av} $ by the Monte Carlo method ** [1] **.

(3) ** Multiple integral calculation by the mean method. **As an example, $\int_1^2 \int_1^2 \frac{1}{x+y} dx =10 ln2-6ln3 \sim 0.33979807 $ To calculate.

Code (1): ** Calculation by simple Monte Carlo method **

simple_Monte_Carlo.py

import numpy as np
from random import random
"""
Simple Monte Carlo integration method:Number of trials N from 10 to 10^Change up to 7
"""
def f(x):
    return 1.0/(1.0+x**2) #Function definition

N_calc_list = [10, 10**2, 10**3, 10**4, 10**5, 10**6, 10**7]

for N in  N_calc_list:
    count = 0.0
    for i in range(N):
        x = random()  # [0,1]Store uniform random numbers up to x in x
        y = random()  # [0,1]Store uniform random numbers up to in y
        if y < f (x):   #If you enter "Mato", count it
            count +=1.0
    area = 4*count/N #Integral result. Evaluation of π

    print(N, ", ", area, ", ", abs((np.pi-area)/np.pi)) #Result output

Result (1)

10 ,  4.0 ,  0.2732395447351627
100 ,  3.08 ,  0.01960555055392467
1000 ,  3.184 ,  0.01349867760918959
10000 ,  3.12 ,  0.006873155106573032
100000 ,  3.14716 ,  0.0017721414021786754
1000000 ,  3.143488 ,  0.0006033075001118286
10000000 ,  3.1414272 ,  5.266551333578959e-05

From the left, the total number of random numbers used in the Monte Carlo method N, integral value, relative error (error from π)

Code (2): ** Monte Carlo integration calculation by mean method **

The_mean_value_method_Monte_Carlo.py

import numpy as np
from numpy.random import rand
"""
Monte Carlo integration by the mean method
The mean value method
"""

    
N_calc_list = [10, 10**2, 10**3, 10**4, 10**5, 10**6, 10**7]
for N in  N_calc_list:
    x_array=rand(N)
    y_array=4.0/(1.0+x_array**2)
    area = (1.0-0.0)*np.sum(y_array)/(N) # y_As av, Σ_i yi/Calculate N and integrate it with the interval(1-0=1)Is hanging

    print(N, ", ", area, ", ",  abs((np.pi-area)/np.pi))#Result output

Result (2)

10 ,  3.18542467193 ,  0.0139521647705
100 ,  3.03821388912 ,  0.0329064827523
1000 ,  3.14697964989 ,  0.0017147341794
10000 ,  3.14560900784 ,  0.00127844526391
100000 ,  3.14380423975 ,  0.000703969738006
1000000 ,  3.14195518509 ,  0.000115397359081
10000000 ,  3.14140220827 ,  6.06206294417e-05

From the left, the total number of random numbers used in the Monte Carlo method N, integral value, relative error (error from π)

Code (3): Double integral by mean method

double_integral.py

import numpy as np
from numpy.random import rand
"""
Monte Carlo integration by the mean method(Double integral)
"""


N_calc_list = [10, 10**2, 10**3, 10**4, 10**5, 10**6, 10**7]
for N in  N_calc_list:
    x_array=rand(N)+1.0 #section[1,2]Uniform random number sequence of x_Store in array
    y_array=rand(N)+1.0 #section[1,2]Uniform random number sequence of y_Store in array
    z_array=1.0/(x_array+y_array)  # z=1/(x+y)Column of z_Store in array
    area = (2.0-1.0)*(2.0-1.0)*np.sum(z_array)/(N)  #Integral calculation
    print(N, ", ", area) #Result output

Result (3)

10 ,  0.346923895691
100 ,  0.343570822229
1000 ,  0.339161537421
10000 ,  0.340241114706
100000 ,  0.339920052165
1000000 ,  0.339757274305
10000000 ,  0.33977097358

From the left, the total number of random numbers used in the Monte Carlo method N, the integral value. The exact solution is $ 10 ln2-6ln3 \ sim 0.33979807 $.

References

[1] Mark Newman, "Computational Physics",Chap10,CreatespaceIndependentPublishingPlatform(2012)

A comment to non-Japanese persons: The two codes above described are both simple for performing numerical integrations using the Monte Carlo methods. One can easily understand how the codes work. The detailed information is accessible via ref. [1].