FEUILLE DE TP 4

---

Intégration numérique (Méthodes composées)

---

Dans ce TP, nous programmons et nous étudions les méthodes composées d'intégration numérique basées sur la méthode d'interpolation de Lagrange.

%matplotlib inline
%config InlineBackend.figure_format = 'retina'

import numpy as np
from scipy.integrate import quad
import matplotlib.pyplot as plt

color_bleu = 'xkcd:dusky blue'
color_vert = 'xkcd:grassy green'
color_orange = 'xkcd:orangish'
color_rouge = 'xkcd:blood orange'
color_bleuclair = 'xkcd:sky blue'
color_vertclair = 'xkcd:grass green'
color_marron = 'xkcd:burnt sienna'
color_violet = 'xkcd:tangerine'

liste_color = [
    color_bleu,
    color_vert,
    color_orange,
    color_rouge,
    color_bleuclair,
    color_vertclair,
    color_marron,
    color_violet
]

Dans ce TP, nous allons programmer de différentes méthodes composées d'intégration numérique. Nous les testerons sur les fonctions suivantes : \begin{align} f_1(x) &= \exp(-x^2),& f_2(x) &= \sqrt{\sin x},\\ %f_2(x) &= 0.1 \,(1+x/3)^3 \log(1+x),& %f_3(x) &= \sin(2\pix) / (1.1 - \sin(\pix)). \end{align}

Question 1.

Définissez ces fonctions et tracez les sur l'intervalle $[0, 3]$.

def f1(x): 
    return np.exp(-x**2)

def f2(x): 
    return np.sqrt(np.sin(x))

liste_f = [f1, f2]
a, b = 0, 3

x = np.linspace(a, b, 200)
fig = plt.figure(figsize=(4*len(liste_f), 4))
fig.suptitle("Les fonctions que l'on va utiliser", fontsize=16)
for k, f in enumerate(liste_f):
    ax = fig.add_subplot(1, len(liste_f), k+1)
    ax.plot(x, f(x), color=color_orange, linewidth=2)
    ax.set_title(f"${f.__name__}$", color=color_marron)

Etant donnée un intervalle $ [ a, b ] $, on considère une subdivision de cet intervalle en $ m $ sous-intervalles de longueur $ h = \frac{ b - a }{ m } $: $ a = a_0 < a_1 < a_2 < \ldots < a_m = b $ avec $ a_j = a + j h $, $ j = 0, 1, \ldots, m $.

Méthode des rectangles à gauche

La formule des rectangles à gauche consiste à faire l'approximation

\begin{equation} I = \int_{a}^{b} f( x ) dx \approx h \sum_{ j = 0 }^{ m - 1 } f( a_i ). \end{equation}

Question 2.

Sur une feuille de papier, tracer la fonction $f_1$ sur l'intervalle $[0, 3]$ et illustrer la méthode des rectangles à gauche pour la subdivision $0 < 0.5 < 1 < 1.5 < 2 < 2.5 < 3$. En déduire une approximation de $\int_{0}^{3} f_1$.

(Vous pouvez également faire ce dessin avec Python à l'aide de la fonction bar de matplotlib.pyplot.)

m = 6
subdivision, pas = np.linspace(a, b, m, endpoint=False, retstep=True)
x = np.linspace(a, b, 1000)

fig = plt.figure(figsize=(4*len(liste_f), 4))
fig.suptitle('Illustration de la méthode des rectangles à gauche')

for k, f in enumerate(liste_f):
    ax = fig.add_subplot(1, len(liste_f), k+1)
    ax.plot(x, f(x), color=color_bleu, linewidth=2)
    ax.bar(
        subdivision, f(subdivision),
        align='edge', width=pas, alpha=0.25, edgecolor=color_bleu
    )
    ax.set_title(f"${f.__name__}$", color=color_bleu)
    
fig = plt.figure(figsize=(4*len(liste_f), 4))
fig.suptitle('Illustration de la méthode des rectangles à droite')

for k, f in enumerate(liste_f):
    ax = fig.add_subplot(1, len(liste_f), k+1)
    ax.plot(x, f(x), color=color_bleu, linewidth=2)
    ax.bar(
        subdivision+pas, f(subdivision+pas),
        align='edge', width=-pas, alpha=0.25, edgecolor=color_bleu
    )
    ax.set_title(f"${f.__name__}$", color=color_bleu)

Approximation : $\int_{0}^{3} f_1(x) dx \approx \frac{1 + 0.8 + 0.4 + 0.1}{2} \approx 1.15$

Iap = 0.
for i in range(6):
    Iap += f1(.5*i)
Iap *= .5
print(Iap)

1.1361627709148368

Question 3.

1. En utilisant la fonction sum du module numpy, programmez une fonction quad_rectg(f,x) qui prend en argument une fonction f et un ndarray x qui contient une subdivision $a_0 < \ldots <a_m$ d'un intervalle $[a,b]$, et retourne la valeur approchée de $I$ obtenue par la formule des rectangles à gauche.
2. Tester ensuite votre méthode sur $f_1$ avec la subdivision de la question 2, et vérifier que le résultat est cohérent avec votre approximation grossière.

def quad_rectg(f, x):
    """
    méthode des rectangles à gauche
    
    Parameters
    ----------
    
    f: function
        la fonction que l'on doit intégrer
    x: ndarray
        la discrétisation choisie
        
    Returns
    -------
    
    Ia: float
        la valeur approchée de l'intégrale
        
    Remarks
    -------
    
    Ia = \sum_{i=1}^{x.size-1} (x_i-x_{i-1}) f(x_{i-1})
    """
    xp, xm = x[1:], x[:-1]
    return np.sum((xp-xm)*f(xm))

m = 6
subdivision = np.linspace(a, b, m+1)
print(quad_rectg(f1, subdivision))

1.1361627709148368

taille = 40
for f in liste_f:
    print("-"*taille)
    print(f"{f.__name__:^{taille}}")
    print(f"Intégrale exacte    : {quad(f, a, b, epsabs=1.e-14)[0]:10.7f}")
    N = 100
    x = np.linspace(a, b, N+1)
    print(f"Rectangles à gauche : {quad_rectg(f, x):10.7f}")
    print("-"*taille)

----------------------------------------
                   f1                   
Intégrale exacte    :  0.8862073
Rectangles à gauche :  0.9012054
----------------------------------------
----------------------------------------
                   f2                   
Intégrale exacte    :  2.3607862
Rectangles à gauche :  2.3539723
----------------------------------------

Question 4.

1. Comparez la valeur approchée de $I$ obtenue à l'aide de la fonction quad_rectg avec sa valeur calculée à l'aide de la fonction quad de scipy.integrate (en choissisant une erreur absolue epsabs = 1e-14) pour des subdivisions de différents pas $h = 1/2, 1/4,1/8, \ldots, 2^{-8}$. Que se passe-t-il quand $h$ diminue ?
2. Tracer l'erreur entre l'intégrale approchée et la valeur "exacte" (obtenue via quad) en fonction de $h$ en échelle log-log. Qu'observez-vous ? Calculer la pente de la droite obtenue.

def analyse(liste_methodes):
    """analyse des méthodes d'intégration"""
    taille = 5
    fig = plt.figure(figsize=(taille*len(liste_f), taille))
    vh = np.array([2**k for k in range(-1, -12, -1)])

    for numf, f in enumerate(liste_f):
        Iex = quad(f, a, b, epsabs=1.e-14)[0]
        ax = fig.add_subplot(1, len(liste_f), numf+1)
        ax.set_xscale('log', base=2)
        ax.set_yscale('log', base=2)

        for k, methode in enumerate(liste_methodes):
            Iap = np.empty(vh.shape)
            for i, h in enumerate(vh):
                x = np.arange(a, b+.5*h, h)
                Iap[i] = methode(f, x)
            ve = abs(Iex - Iap)
            ax.scatter(
                vh, ve, c=liste_color[k], s=50,
                label=f"{methode.__name__}"
            )
            indices = ve > 1.e-14
            p = np.polyfit(
                np.log(vh[indices]), np.log(ve[indices]), 1
            )[0]
            ax.plot(vh, ve[0]*(vh/vh[0])**p, label=f'pente {p:5.3f}', color=liste_color[k], linestyle=(0, (p, 1)))

        ax.set_title(f"Fonction $f_{numf}$", fontsize=20)
        ax.set_xlabel(r'$h$', fontsize=16)
        ax.set_ylabel(r'$e_h$', fontsize=16)
        ax.legend()

    fig.tight_layout()

analyse([quad_rectg])

Méthode des trapèzes

On rappelle que la formule des trapèzes consiste à faire l'approximation

\begin{equation} \int_{a}^{b} f( x ) dx \approx h \sum_{j=0}^{m-1} \frac{1}{2} \Bigl( f(a_j)+f(a_{j+1}) \Bigr) =h \Bigl[\frac{1}{2} \Bigl( f(a)+f(b) \Bigr) + \sum_{j=1}^{m-1} f(a_j) \Bigr] \end{equation}

Question 5.

1. Reprenez la question 2 avec la méthode des trapèzes.
2. Programmez une fonction quad_trapeze(f,x) qui prends un argument une fonction f et un ndarray x qui contient une subdivision $a_0 < \ldots <a_m$ d'un intervalle $[a,b]$, et retourne la valeur approchée de $I$ obtenue par la formule des trapèzes.
3. Comparez la méthode des trapèzes à la méthode des rectangles en traçant sur une même figure les erreurs en fonction du pas $h$ en echelle log-log. Laquelle est la plus satisfaisante ?

def quad_trapeze(f, x):
    """
    méthode des trapèzes
    
    Parameters
    ----------
    
    f: function
        la fonction que l'on doit intégrer
    x: ndarray
        la discrétisation choisie
        
    Returns
    -------
    
    Ia: float
        la valeur approchée de l'intégrale
        
    Remarks
    -------
    
    Ia = \sum_{i=1}^{x.size-1} .5*(x_i-x_{i-1}) (f(x_{i-1})+f(x_i))
    """
    xp, xm = x[1:], x[:-1]
    return .5*np.sum((xp-xm)*(f(xm)+f(xp)))

analyse([quad_rectg, quad_trapeze])

Méthode du point milieu

On rappelle que la formule du point milieu consiste à faire l'approximation

\begin{equation} \int_{a}^{b} f(x)dx \approx h \sum_{j=0}^{m-1} f(a_{j+{1}/{2}}), \quad \text{où} \quad a_{ j+{1}/{2} } = \frac{ a_j + a_{j+1} }{ 2 }. \end{equation}

Question 6.

1. Programmez une fonction quad_milieu(f,x) qui prends un argument une fonction f et un ndarray x qui contient une subdivision $a_0 < \ldots <a_m$ d'un intervalle $[a,b]$, et retourne la valeur approchée de $I$ obtenue par la formule du point milieu.
2. Comparez cette méthode à celle des trapèzes en traçant sur une même figure les erreurs en fonction du pas $h$ en echelle log-log. L'une des deux méthodes est-elle plus intéressante ?

def quad_milieu(f, x):
    """
    méthode des points milieux
    
    Parameters
    ----------
    
    f: function
        la fonction que l'on doit intégrer
    x: ndarray
        la discrétisation choisie
        
    Returns
    -------
    
    Ia: float
        la valeur approchée de l'intégrale
        
    Remarks
    -------
    
    Ia = \sum_{i=1}^{x.size-1} (x_i-x_{i-1}) f(.5*(x_{i-1}+x_i))
    """
    xp, xm = x[1:], x[:-1]
    return np.sum((xp-xm)*f(.5*(xm+xp)))

analyse([quad_rectg, quad_trapeze, quad_milieu])

Méthode de Simpson

On rappelle la formule de Simpson qui consiste à faire l'approximation

\begin{equation} \int_{a}^{b} f( x ) dx \approx h \sum_{j=0}^{m-1} \frac{1}{6} \Bigl( f(a_j)+4f(a_{j+{1}/{2}})+f(a_{j+1}) \Bigr). \end{equation}

Question 7.

1. Programmez une fonction quad_simpson(f,x) qui prend en argument une fonction f et un ndarray x qui contient une subdivision $ a_0 < \ldots < a_m $ d'un intervalle $[a,b]$, et retourne la valeur approchée de $I$ obtenue par la formule de Simpson. Comparer avec les méthodes précédentes et commenter.

def quad_simpson(f, x):
    """
    méthode de Simpson
    
    Parameters
    ----------
    
    f: function
        la fonction que l'on doit intégrer
    x: ndarray
        la discrétisation choisie
        
    Returns
    -------
    
    Ia: float
        la valeur approchée de l'intégrale
        
    Remarks
    -------
    
    Ia = \sum_{i=1}^{x.size-1} (x_i-x_{i-1})/6 * (
        f(x_{i-1}) + 4 f(.5*(x_{i-1}+x_i)) + f(x_i)
    )
    """
    return (quad_trapeze(f, x) + 2*quad_milieu(f, x))/3

analyse([quad_rectg, quad_milieu, quad_trapeze, quad_simpson])