Intégration
===========

Soit une fonction :math:`f:\mathbb{R}\rightarrow\mathbb{R}` dont on cherche à calculer l'intégrale suivante :

.. math::
    I = \int_a^bf(x)\,dx
    
Le calcul numérique de cette intégrale est nécessaire dans les cas suivants :

1. L'expression de :math:`f` est connue mais sa primitive ne l'est pas.
2. Le calcul de :math:`f(x)` se fait par une simulation et la primivive n'est pas connue.
3. Seules les valeurs de :math:`f` pour un échantillonnage d'un intervalle :math:`[a,b]` sont connues, par exemple des valeurs expérimentales.

Méthode des rectangles
----------------------

Considérons l'échantillonnage suivant de la fonction sur l'intervalle :math:`[a,b]`, pour :math:`n=0,1,\ldots N-1` :

.. math::
    \begin{align}
    &x_n = a+\frac{b-a}{N-1}n\\
    &y_n = f(x_n)
    \end{align}
    
La méthode des rectangles consiste à utiliser l'approximation suivante :

.. math::
    \int_{x_n}^{x_{n+1}}f(x)dx\approx f(x_n)(x_{n+1}-x_n)
    
La figure suivante en montre le principe graphiquement : l'intégrale sur chaque sous-intervalle est assimilée à l'aire d'un rectangle, c'est-à-dire que la fonction est assimilée à une constante égale à sa valeur à gauche.
    
.. image:: methodeRectangle.svg
    :scale: 80
    
Cette méthode est nommée **méthode des rectangles à gauche** car, pour chaque sous-intervalle, on utilise la valeur à gauche du sous-intervalle. L'approximation de l'intégrale s'écrit :

.. math::
    \int_a^bf(x)\,dx\approx \sum_{n=0}^{N-2}f(x_n)(x_{n+1}-x_n)=\frac{b-a}{N-1}\sum_{n=0}^{N-2}f(x_n)
    :label: rectangles
    
La limite de cette expression lorsque :math:`N\rightarrow\infty` est précisément la définition de l'intégrale de Riemann, ce qui garantit la convergence vers l'intégrale lorsque :math:`N` augmente, c'est-à-dire lorsque le pas 

.. math::
    h=\frac{b-a}{N-1}
    
est de plus en plus petit. L'expression :eq:`rectangles` montre que le calcul de l'intégrale par la méthode des rectangles revient en fait à calculer la somme suivante :

.. math::
    S = \sum_{n=0}^{N-2}y_n
    
Pour la méthode des rectangles à droite, on aurait à calculer la somme 

.. math::
    S' = \sum_{n=1}^{N-1}y_n
    
Lorsque :math:`N` est grand, ces deux sommes sont très proches et les deux méthodes sont donc équivalentes.

Intégration d'une fonction
--------------------------

Supposons qu'on dispose d'une expression de la fonction :math:`f` et qu'une fonction Python évalue cette expression. La signature d'une telle fonction est :

.. code-block:: python

    def f(x):
        # calcul et renvoie de f(x)
        
La fonction suivante réalise le calcul approché par la méthode des rectangles (à gauche) :

.. code-block:: python

    import numpy as np

    def integrale(f,a,b,N):
        X = np.linspace(a,b,N)
        Y = f(X)
        return (b-a)/(N-1)*np.sum(Y[0:N-1])
    
Pour calculer la somme, nous sélectionnons les éléments d'indice 0 à N-2 au moyen d'un tranche puis nous appelons la fonction :py:obj:`np.sum` sur le tableau résultant. Il est bien sûr possible d'écrire explicitement la boucle de calcul de la somme :

.. code-block:: python

    def integrale_moins_rapide(f,a,b,N):
        X = np.linspace(a,b,N)
        Y = f(X)
        S = 0
        for n in range(N-1):
            S += Y[n]
        return (b-a)/(N-1)*S
        
Pour tester une méthode numérique, il est judicieux de le faire sur un cas où la solution exacte est connue. Considérons l'intégrale de la fonction sinus sur l'intervalle :math:`[0,\pi/2]`, qui vaut exactement 1 :

.. code-block:: python
    
    def f(x):
        return np.sin(x)

    print(integrale(f,0,np.pi/2,100))
    print(integrale(f,0,np.pi/2,1000))
    print(integrale(f,0,np.pi/2,10000))
    
    0.9920457059690396
    0.9992136096236359
    0.9999214502723134

Ces résultats suggèrent que l'erreur est divisée par 10 lorsque :math:`N` est augmenté d'un facteur 10. La méthode des rectangles est en effet d'ordre 1, c'est-à-dire que la valeur absolue de l'erreur est majorée par :math:`K|h|`, où :math:`K` est une constante.
    
Fonction intégrale
------------------

Il arrive fréquemment qu'on ait besoin de représenter graphiquement une fonction qui s'exprime sous la forme d'une intégrale :

.. math::
    g(x)=\int_0^xf(x')\,dx'
    
Il serait bien sûr extrêmement inefficace d'appeler la fonction :py:obj:`integrale` (définie plus haut) pour chaque valeur de :math:`x` de l'intervalle sur lequel on veut représenter la fonction. La bonne méthode consiste à échantillonner l'intervalle sur lequel on veut calculer :math:`g(x)` puis calculer l'intégrale de proche en proche au moyen de la relation :

.. math::
    \begin{align}
    g(x_{n+1}) = &g(x_{n})+\int_{x_{n}}^{x_{n+1}}f(x')\,dx'\\
               \approx &g(x_{n})+f(x_n)(x_{n+1}-x_n)
    \end{align}

Ce schéma de calcul consitue en fait la **méthode d'Euler** apliquée à l'équation différentielle suivante :

.. math::
    \begin{align}
        &\frac{dg}{dx}=f(x)\\
        &g(0)=0
    \end{align}
    
Voyons par exemple la fonction *erreur*, définie par :

.. math::
    {\rm erf}(x)=\frac{2}{\sqrt{\pi}}\int_0^xe^{-t^2}\,dt
    
Cette fonction existe dans le module `scipy.special`__, qui permet d'évaluer un grand nombre de fonctions spéciales, mais nous souhaitons, à titre d'exemple, obtenir sa représentation graphique par la méthode des rectangles.

__ https://docs.scipy.org/doc/scipy/reference/special.html


.. code-block:: python

    import numpy as np
    import matplotlib.pyplot as plt
    
    N = 1000
    x = np.linspace(0,3,N)
    h = x[1]-x[0]
    erf = np.zeros(N)
    I = 0
    for n in range(N):
        erf[n] = I
        I += np.exp(-x[n]**2)
    erf *= 2*h/np.sqrt(np.pi)
    
    plt.figure()
    plt.plot(x,erf)
    plt.xlabel('x',fontsize=16)
    plt.ylabel('erf(x)',fontsize=16)
    plt.grid()
    plt.show()
    
.. image:: fonctionErreur.png
    :scale: 80