Cours 4 : Monte-Carlo vs rectangles en dimension 5

Comparer le temps de calcul de la méthode des rectangles et de la méthode de Monte Carlo en dimension 5 avec la fonction f donnée.

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import time

np.__version__, matplotlib.__version__

('1.17.2', '3.1.1')

# fonction dont il faut calculer l'integrale
def f(u, v, w, x, y):
    return np.exp(u) * 2*v * np.sin(w) * np.cos(x) * (1 / (y+1))

# valeur de l'integrale 
a, b = 0, 1
integrale_f = (
    (np.exp(b) - np.exp(a)) *
    (b**2 - a**2) *
    (-np.cos(b) + np.cos(a))  *
    (np.sin(b) -  np.sin(a)) * 
    (np.log(1+b) - np.log(1+a))
    )
integrale_f

0.46071391454414323

Methode des rectangles pour le calcul de l'intégrale de f

n_points, a, b = 10, 0, 1

def calcul_methode_rectangle_boucle_for():
    integrale_rec = 0
    
    # on definit la grille 1D puis 5D
    grid_1D = np.linspace(a, b, n_points, endpoint=False)
    step = grid_1D[1] - grid_1D[0]
    grid_1D += 0.5*step 
    u_grid, v_grid, w_grid, x_grid ,y_grid = np.meshgrid(*[grid_1D.copy() for _ in range(5)])
    
    # on "applatit" les grilles pour les parcourir avec une seule boucle for (plutot que 5)
    u_grid, v_grid, w_grid, x_grid ,y_grid = u_grid.flatten(), v_grid.flatten(), w_grid.flatten(), x_grid.flatten() ,y_grid.flatten()

    # on parcours la grille avec une boucle for pour calculer l'integrale
    for i in range(u_grid.size):
        u, v, w, x, y = u_grid[i], v_grid[i], w_grid[i], x_grid[i] ,y_grid[i]
        integrale_rec += f(u, v, w, x, y) * step**5


    return integrale_rec

def calcul_methode_rectangle_vectorielle():
    
    # on definit la grille 1D puis 5D
    grid_1D = np.linspace(a, b, n_points, endpoint=False)
    step = grid_1D[1] - grid_1D[0]
    grid_1D += 0.5*step 
    grid_5D = np.meshgrid(*[grid_1D.copy() for _ in range(5)])
    # calcul direct de l'integrale
    return f(*grid_5D).sum() * step**5

print('Methode 1: boucle for')
%time valeur = calcul_methode_rectangle_boucle_for()
print(f'difference : {valeur - integrale_f}')

print('\nMethode 2: vetorielle')
%time valeur = calcul_methode_rectangle_vectorielle()
print(f'difference : {valeur - integrale_f}')

Methode 1: boucle for
CPU times: user 869 ms, sys: 19 ms, total: 888 ms
Wall time: 1 s
difference : -1.5291653360782753e-05

Methode 2: vetorielle
CPU times: user 7.8 ms, sys: 1.18 ms, total: 8.98 ms
Wall time: 6.89 ms
difference : -1.529165336255911e-05

Methode de Monte-Carlo pour le calcul de l'intégrale de f

# echantilloner variable uniforme sur [a,b]^5 
N = n_points**5

def calcul_methode_mc_boucle_for():
    integrale_mc = 0
    # echantilloner une loi uniforme sur [a,b]^5
    uniform = a + (b-a) * np.random.rand(N, 5)

    # parcourir les valeurs echantillonnees pour calculer la somme de monte carlo
    for i in range(N):
        u, v, w, x, y = uniform[i]
        integrale_mc += f(u, v, w, x, y) / N
        
    return integrale_mc

def calcul_methode_mc_vectorielle():
    # echantilloner 5 lois uniformes sur [a,b] pour chacune des coordonnees
    uniform = [a + (b-a) * np.random.rand(N) for _ in range(5)]

    # calcul direct  
    return f(*uniform).sum() / N

print('Methode 1: boucle for')
%time valeur = calcul_methode_mc_boucle_for()
print(f'difference : {valeur - integrale_f}')

print('\nMethode 2: vectorielle')
%time valeur = calcul_methode_mc_vectorielle()
print(f'difference : {valeur - integrale_f}')

Methode 1: boucle for
CPU times: user 1.14 s, sys: 19.7 ms, total: 1.16 s
Wall time: 1.07 s
difference : 0.0011634917967472402

Methode 2: vectorielle
CPU times: user 11.5 ms, sys: 315 µs, total: 11.9 ms
Wall time: 10.2 ms
difference : 0.003263222379336128