## Méthode des Trapèzes



On sait que l'aire d'un trapèze est donnée par : $$Aire = \frac{(B + b)h}{2}$$

On définit : $$W_i = \frac{(f(x_i) + f(x_{i+1}))h}{2}$$

On trouve une formule pour approximer notre intégrale :

$$
\begin{align*}
    I &\approx \sum_{i=0}^{n-1} W_i \\
      &\approx \sum_{i=0}^{n-1} \frac{h}{2}(f(x_i) + f(x_{i+1}) \\ 
      &\approx \frac{h}{2}(f(x_0) + f(x_1) + f(x_1) + ... + f(x_{n-1}) + f(x_{n-1}) + f(x_n)) \\
      &\approx \frac{h}{2} \left[ f(x_0) + 2 \sum_{i=0}^{n-1}f(x_i) + f(x_n)) \right] \\
      &\approx h \left[ \frac{1}{2}(f(x_0) + f(x_n)) + \sum_{i=0}^{n-1}f(x_i)) \right]
\end{align*}
$$

### Itération 
- On diminue progressivement le pas de h de $\frac{1}{2}$ à chaque itération jusqu'à atteindre la précision souhaitée.
- /!\ une itération = une évaluation complète de l'intégrale.

Pour une fonction : $$\int_a^b f(x)dx$$

$n = 2^{k-1}$ => 1, 2, 4, 8,...

$h = \frac{b-a}{n}$ => $h$, $\frac{h}{2}$, $\frac{h}{4}$,...

Inconvénient : convergence lente donc besoin de beaucoup d'itérations pour être précis.


```python
import numpy as np
import math

# precision de l'évaluation numérique de l'intégrale (trapèzes)
TRESHOLD_PRECISION = 1e-9
DIMENSIONS = [1] #[1, 2, 10, 15]
X_0 = -3 
X_N = 3

class Trapezes():
  def __init__(self, f, x_0, x_n):
    self.f = f
    self.x_0 = x_0
    self.x_n = x_n
  
  # borne inférieure : x_0
  # borne supérieure : x_n
  # en pratique on itère sur la méthode 
  # et on diminue le pas h progressivement (de 1/2 à chaque itération)
  # k : itération en cours
  def trapezes(self, k):
    # si k = 1, le nombre de trapèzes serait de 1 (en effet, 2^0 = 1)
    if k == 1: 
      # Aire d'un trapèze
      h = (self.x_n - self.x_0)
      return (self.f(self.x_0) + self.f(self.x_n)) * h / 2

    # nombre de trapèzes
    n = int(math.pow(2, k-1))
    # intervalle h
    h = (self.x_n - self.x_0) / n

    sum = 1/2 * (self.f(self.x_0) + self.f(self.x_n))
    for i in range(1, n):
      x_i = self.x_0 + (i * h)
      sum += self.f(x_i)

    return h * sum 
  
  # N : nombre de dimensions
  def integrate(self, iterations, precision):
    new_integral = 0
    for k in range(1, iterations):
      old_integral = new_integral
      new_integral = self.trapezes(k)
      if abs(old_integral - new_integral) < precision and k > 1:
        return new_integral, k
        break

def gaussian(x):
  f = math.exp(-(x**2) / 2)
  return f

def gaussian_2_dim(x, y):
  return math.exp(-(x**2 + y**2) / 2)

def gaussian_theo(N):
  return math.sqrt((2*math.pi)**N)

print(f"dimension 1")
print("-----------")
print(f"Valeur théorique approximative: {gaussian_theo(N=1)}")

gaussian_trapezes = Trapezes(gaussian, X_0, X_N)
gaussian_trapezes_value = gaussian_trapezes.integrate(iterations = 20, precision = TRESHOLD_PRECISION)
print(f"Valeur de la gaussienne à 1 dimension (trapèzes): {gaussian_trapezes_value}")

print(f"dimension 2")
print("-------------")
print(f"Valeur théorique approximative: {gaussian_theo(N=2)}")

print("-------------")

def universe_model(x, omega_m, omega_rad, omega_k, omega_lambda):
    

gaussian_trapezes = Trapezes(gaussian, X_0, X_N)
gaussian_trapezes_value = gaussian_trapezes.integrate(iterations = 20, precision = TRESHOLD_PRECISION)
print(f"Valeur de la gaussienne à 1 dimension (trapèzes): {gaussian_trapezes_value}")

```

    dimension 1
    -----------
    Valeur théorique approximative: 2.5066282746310002
    Valeur de la gaussienne à 1 dimension (trapèzes): (2.4998608892968623, 16)
    dimension 2
    -------------
    Valeur théorique approximative: 6.283185307179586


## Méthode de Romberg

Combine la méthode itérative des trapèzes avec l'extrapolation de Richardson.    
**Erreur sur la méthode des Trapèzes**: $E(h) = c_1h^2 + c_2h^4 + ... = \sum_i c_ih^{2i}$

### Extrapolation de Richardson

Pour un pas de h : 
$$G = I(h) + E(h)$$

Pour un pas de h/2 : 
$$G = I(h/2) + E(h/2)$$

On obtient le système :

$$
\begin{align}
\begin{cases}
    G = I(h) + ch^p \\ 
    G = I(h/2) + c(h/2)^p \implies 2^pG = 2^pI(h/2) + ch^p
\end{cases}
\end{align}
$$

On soustrait l'équation $(1)$ à l'équation $(2)$ :

$$
\begin{align}
    2^pG - G &= 2^pI(h/2) - I(h) \\
    (2^p - 1)G &= 2^p I(h/2) - I(h) \\
    G &= \frac{2^p I(h/2) - I(h)}{2^p - 1}
\end{align}
$$

Nuos avons donc une formule pour l'extrapolation de Richardson : 

$$
G = \frac{2^p I(h/2) - I(h)}{2^p - 1}
$$

**Concrètement**:

On va se constituer une matrice $n x n$ avec _n : le nombre d'itérations_. Elle contiendra dans sa première colonne les approximations de l'intégrale par la méthode des trapèzes. Ensuite on rempli "le reste" en cascade avec la formule de l'extrapolation de Richardson. 

**Algorithme**:

On remarque que plus le nombre d'itérations augmente, plus la précision augmente.




```python
import numpy as np
import math

class Romberg():
  def __init__(self, f, x_0, x_n):
    self.f = f
    self.x_0 = x_0
    self.x_n = x_n

  def trapezes(self, k):
    # si k = 1, le nombre de trapèzes serait de 1 (en effet, 2^0 = 1)
    if k == 1: 
      # Aire d'un trapèze
      h = (self.x_n - self.x_0)
      return (self.f(self.x_0) + self.f(self.x_n)) * h / 2

    # nombre de trapèzes
    n = int(math.pow(2, k-1))
    # intervalle h
    h = (self.x_n - self.x_0) / n

    sum = 1/2 * (self.f(self.x_0) + self.f(self.x_n))
    for i in range(1, n):
      x_i = self.x_0 + (i * h)
      sum += self.f(x_i)

    return h * sum 

  def richardson_extrapolation(self, I_prev, I_current, k):
    return ((2**k * I_current) - I_prev) / (2**k - 1)

  def integrate(self, iterations, precision):
    # R(n, m)
    self.map = np.zeros((iterations, iterations))

    """
    iterations = 3
       j
       __________________
    i | trap   x    x    |
      | trap rich   x    |
      | trap rich approx |  
       ------------------
    """
    for i in range(1, iterations):
      i_prime = i - 1
      # Méthode des trapèzes
      self.map[i_prime, 0] = self.trapezes(i)
      for j in range(1, i):
        # Extrapolation de Richardson
        self.map[i_prime, j] = self.richardson_extrapolation(self.map[i_prime-1, j-1], self.map[i_prime, j-1], 2*j)
      
      if abs(self.map[i_prime - 1, i_prime - 1] - self.map[i_prime, i_prime]) < precision:
        return self.map[i_prime, i_prime], i
        break

# precision de l'évaluation numérique de l'intégrale (trapèzes)
TRESHOLD_PRECISION = 1e-9
DIMENSIONS = [1] #[1, 2, 10, 15]
X_0 = -3 
X_N = 3

def gaussian(x):
  f = math.exp(-(x**2) / 2)
  return f
        
# romberg
gaussian_romberg = Romberg(gaussian, X_0, X_N)
gaussian_romberg_value = gaussian_romberg.integrate(iterations = 10, precision = TRESHOLD_PRECISION)
print(f"Valeur de la gaussienne à 2 dimension (romberg): {gaussian_romberg_value}")
```

    Valeur de la gaussienne à 2 dimension (romberg): (2.4998608894818006, 8)


### Méthode de Monte-Carlo

- Plus efficace pour calculer des intégrales à plusieurs dimensions.
- On tire les points aléatoirement dans le domaine de la fonction.

$$
\begin{align}
    <f> &= \frac{1}{b-a} \int_a^b f(x)dx \\
    \int_a^b f(x)dx &= (b-a)<f> \\
    \int_a^b f(x)dx &\approx (b-a) \frac{1}{N} \sum_{i=0}^{n-1} f(x_i)
\end{align}
$$


```python
import math
import random
import numpy as np

def gaussian_theo(N):
  return math.sqrt((2*math.pi)**N)

def gaussian(x_vec):
    sum = 0
    # x^2 + y^2 + z^2 + ...
    for x_i in x_vec:
        sum += x_i**2
        
    return math.exp(-(sum) / 2)

class MonteCarlo():
    def __init__(self, f, x_0, x_n):
        self.f = f
        self.x_0 = x_0
        self.x_n = x_n
    
    def integrate(self, iterations):
        sum = 0
        # e.g.: 3 = intégrale triple
        dimension = len(self.x_0) 
        # N * N * N * ... = N^(dimension)
        volume = np.prod([x_n - x_0 for x_0, x_n in zip(self.x_0, self.x_n)]) 
        # somme des f(x_i)
        sum = 0
        for i in range(iterations):
            x_i = [random.uniform(self.x_0[k], self.x_n[k]) for k in range(dimension)]
            sum += self.f(x_i)
        
        integral = volume * (1/iterations) * sum
        return integral
    

print(f"predicted gaussian integral value: {gaussian_theo(3)}")
integral = MonteCarlo(gaussian, [-3, -3, -3], [3, 3, 3]).integrate(1000000)
print(f"approximated gaussian integral value: {integral}")
```

    predicted gaussian integral value: 15.749609945722419
    approximated gaussian integral value: 15.627813816362522


#### Amélioration 

- Utiliser une répartition arbitraire selon la forme de la fonction à intégrer.

c2dec1c9faff3285568bb150531554f51a9ae7e5
ipynb
Jupyter Notebook
python/notebooks/integration.ipynb
Mathieu-R/computational-physics
7b1a39ffb65f864bc741a3e709892ddb18723a58
MIT
python/notebooks/integration.ipynb
Mathieu-R/computational-physics
7b1a39ffb65f864bc741a3e709892ddb18723a58
MIT
2021-06-08T23:08:56.000Z
2022-03-12T00:54:45.000Z
python/notebooks/integration.ipynb
Mathieu-R/computational-physics
7b1a39ffb65f864bc741a3e709892ddb18723a58
MIT
Qwen/Qwen-72B
1. YES
2. YES