Calcul de moments d'inertie par intégration de Monte-Carlo

3.a. Principe

La méthode consiste à tirer aléatoirement N points dans le domaine carré Ω avec une densité de probabilité uniforme, égale à l'inverse de l'aire du carré :

$$p(x,y)=\frac{1}{4}\text{(4)}$$

Pour chacun de ces N points, on calcule :

$$f_i=f(x_i,y_i)=(x_i^2+x_i^2)\rho (x_i,y_i)\text{(5)}$$

Cette fonction est très facile à évaluer. Il suffit de faire des tests pour savoir si le point se trouve dans la matière, et de renvoyer x²+y² si c'est le cas, 0 dans le cas contraire.

Une estimation de l'intégrale est obtenue en calculant la moyenne empirique de f/p, c'est-à-dire :

$$S_N=moy(f/p,N)=\frac{1}{N}\sum _{i=0}^{N-1}\frac{f_i}{p_i}\text{(6)}$$

Dans le cas présent, la densité est uniforme (p_i=p), donc on peut diviser par p après le calcul de la somme. On calcule aussi la variance empirique de la manière suivante :

$$v_N=var(f/p,N)=\frac{1}{N}\sum _{i=0}^{N-1}(\frac{f_i}{p_i}{)}^2-S_N^2\text{(7)}$$

L'intervalle de confiance à 95 pour cent pour l'intégrale est :

$$[S_N-\frac{1,96\sqrt{v_N}}{\sqrt{N}},S_N+\frac{1,96\sqrt{v_N}}{\sqrt{N}}]\text{(8)}$$

3.b. Réalisation

import numpy
import numpy.random
from matplotlib.pyplot import *
import math
                

On commence par écrire la fonction qui calcule f(x,y) :

def f_Jz(x,y,param):
    R1 = param[0]
    R2 = param[1]
    a = param[2]
    # ....
                

L'argument param est une liste qui contient les paramètres de la roue. Cette fonction doit être écrite avec le souci de minimiser le nombre de tests et de calculs.

def f_Jz(x,y,param):
    R1 = param[0]
    R2 = param[1]
    a = param[2]
    r2 = x*x+y*y
    if r2 > 1 or r2 < R1*R1:
        return 0.0
    if r2 > R2*R2:
        return r2
    if abs(x) < a:
        return r2
    if abs(y) < a:
        return r2
    return 0.0
                

La fonction de calcul de l'intégrale est :

integration2d(fonction,N,param,np)

Cette fonction fait le calcul de la somme, de la variance et de l'intervalle de confiance à 95 pour cent. Elle génère aussi une image carrée dont le nombre de pixels sur chaque dimension est np. Cette image, noire au départ, est complétée par un pixel blanc pour chaque point tiré dont la valeur de f est non nulle.

def integration2d(fonction,N,param,np):
    image = numpy.zeros((np,np))
    x = -1+2*numpy.random.random_sample(N)
    y = -1+2*numpy.random.random_sample(N)
    p = 1.0/4
    somme = 0.0
    somme2 = 0.0
    for i in range(N):
        f = fonction(x[i],y[i],param)
        somme += f
        somme2 += f*f
        if f!=0:
            image[int((y[i]+1)/2*np),int((x[i]+1)/2*np)] = 1.0
    moyenne = somme/(p*N)
    variance = somme2/(p*p*N)-moyenne*moyenne
    return (moyenne,math.sqrt(variance/N)*1.96,image)
                 

Voici un calcul :

param = [0.1,0.8,0.15]
N = 10**5
(integrale,intervalle,image) = integration2d(f_Jz,N,param,500)
figure()
imshow(image,cmap='gray')
                 

figA.pdf

print((integrale,intervalle))
--> (1.1361162182318325, 0.009031678004722264)

4.a. Échantillonnage sur le disque

Une première méthode consiste à échantillonner uniformément sur le disque de rayon 1. On doit pour cela utiliser les coordonnées polaire (r,θ) et faire des tirages avec la densité de probabilité uniforme sur le disque :

$$p(r,\theta )=\frac{1}{\pi }\text{(9)}$$

Pour réaliser cet échantillonnage à partir de variables à densité uniforme, on considère la fonction de répartition bidimensionnelle :

$$F(r,\theta )={\int }_0^{\theta }{\int }_0^r\frac{r\text{'}d\theta \text{'}dr\text{'}}{\pi }=\frac{\theta }{2\pi }{r}^2\text{(10)}$$

L'inversion de la fonction de répartition permet d'échantillonner r et θ à partir de deux variables (u₁,u₂) uniformes sur l'intervalle [0,1] :

$$\begin{array}{ccc}\hfill & \theta =2\pi u_1\hfill & \hfill \text{(11)}\\ \hfill & r=\sqrt{u_2}\hfill & \hfill \text{(12)}\end{array}$$

4.b. Réalisation

On redéfinit la fonction f en coordonnées polaires :

def f_Jz(r,theta,param):
    R1 = param[0]
    R2 = param[1]
    a = param[2]
    # .....
                

def f_Jz_polaire(r,theta,param):
    R1 = param[0]
    R2 = param[1]
    a = param[2]
    r2 = r*r
    if r2 > 1 or r2 < R1*R1:
        return 0.0
    if r2 > R2*R2:
        return r2
    if abs(r*math.cos(theta)) < a:
        return r2
    if abs(r*math.sin(theta)) < a:
        return r2
    return 0.0
                

et la fonction d'intégration :

integration2d_disque(fonction,N,param,np)
                

def integration2d_disque(fonction,N,param,np):
    image = numpy.zeros((np,np))
    theta = 2*numpy.pi*numpy.random.random_sample(N)
    r = numpy.power(numpy.random.random_sample(N),0.5)
    somme = 0.0
    somme2 = 0.0
    p = 1.0/numpy.pi
    for i in range(N):
        f = fonction(r[i],theta[i],param)/p
        somme += f
        somme2 += f*f
        if f!=0:
            image[int((r[i]*math.cos(theta[i])+1)/2*np),int((r[i]*math.sin(theta[i])+1)/2*np)] = 1.0
    moyenne = somme/N
    variance = somme2/N-moyenne*moyenne
    return (moyenne,math.sqrt(variance/N)*1.96,image)
                

Voici un calcul :

(integrale,intervalle,image) = integration2d_disque(f_Jz_polaire,N,param,500)
figure()
imshow(image,cmap='gray')
                 

figB.pdf

print((integrale,intervalle))
--> (1.1338844321238779, 0.007306451381173197)

L'écart est bien réduit par rapport à l'échantillonnage sur le domaine carré.

4.c. Échantillonnage non uniforme

Pour réduire encore la variance, on peut chercher à échantillonner sur le disque de rayon 1, avec une densité proportionnelle à x²+y², c'est-à-dire :

$$p(r,\theta )=\frac{2r^2}{\pi }\text{(13)}$$

La fonction de répartition est :

$$F(r,\theta )=\frac{\theta }{2\pi }{r}^4\text{(14)}$$

On en déduit l'échantillonnage à partir de deux variables uniformes sur l'intervalle [0,1] :

$$\begin{array}{ccc}\hfill & \theta =2\pi u_1\hfill & \hfill \text{(15)}\\ \hfill & r=u_2^{\frac{1}{4}}\hfill & \hfill \text{(16)}\end{array}$$

4.d. Réalisation

On écrit la fonction d'intégration :

integration2d_disque2(fonction,N,param,np)
                

def integration2d_disque2(fonction,N,param,np):
    image = numpy.zeros((np,np))
    theta = 2*numpy.pi*numpy.random.random_sample(N)
    r = numpy.power(numpy.random.random_sample(N),0.25)
    somme = 0.0
    somme2 = 0.0
    k = numpy.pi/2
    for i in range(N):
        f = fonction(r[i],theta[i],param)*k/(r[i]*r[i])
        somme += f
        somme2 += f*f
        if f!=0:
            image[int((r[i]*math.cos(theta[i])+1)/2*np),int((r[i]*math.sin(theta[i])+1)/2*np)] = 1.0
    moyenne = somme/N
    variance = somme2/N-moyenne*moyenne
    return (moyenne,math.sqrt(variance/N)*1.96,image)
                

Voici un calcul :

(integrale,intervalle,image) = integration2d_disque2(f_Jz_polaire,N,param,500)
figure()
imshow(image,cmap='gray')
                 

figC.pdf

print((integrale,intervalle))
--> (1.1370994609655842, 0.004352596815941401)

Cet échantillonnage préférentiel permet de réduire la variance de manière importante. Il peut être utilisé pour différentes formes d'objets.

Il est intéressant de comparer avec le moment d'inertie d'une roue pleine (avec seulement le trou central) :

param = [0.1,0.1,0]
(integrale,intervalle,image) = integration2d_disque2(f_Jz_polaire,N,param,500)
figure()
imshow(image,cmap='gray')
                  

figD.pdf

print((integrale,intervalle))
--> (1.5706549551235363, 9.235866842716063e-05)

Dans ce cas la variance est très faible. En effet, un disque plein complet doit donner une variance théorique nulle avec cet échantillonnage. La contribution du trou central au moment d'inertie est très faible.

4.e. Distribution des intégrales évaluées

Obtenir un histogramme donnant la distribution des évaluations pour un ensemble de calculs faits avec un nombre N de tirages donné. Il faudra choisir judicieusement l'intervalle de valeurs et le nombre de sous-intervalles pour le calcul de cet histogramme. Pour des raisons de temps de calcul, la valeur de N sera moins grande que précédemment.