TangenteX.com

L'intégration numérique

L'intégration numérique est un chapitre important de l'analyse numérique et un outil indispensable en physique numérique. On intègre numériquement dans deux cas principaux :

  • on ne peut pas intégrer analytiquement,
  • l'intégrande est fourni non pas sous la forme d'une fonction mais de tableaux de mesures, cas d'ailleurs le plus fréquent dans la vraie vie.

Les méthodes numériques d'intégration d'une fonction sont nombreuses et les techniques très diverses. Des très simples, comme la méthode des rectangles aux très complexes comme certaines variétés de la méthode de Monte-Carlo. Nous n'aborderons ici que des méthodes (ou schémas) simples et courants. Mon but est de vous donner un outil pour intégrer des fonctions pas très tourmentées, celles que l'on rencontre en physique dans le premier cycle universitaire. Pour les autres schémas, voir par exemple la méthode de Monte-Carlo...

La méthode des rectangles

L'algorithme

Considérons donc une fonction de \( \mathbb{R} \) dans \( \mathbb{R} \) continue sur un intervalle [a,b]. Pour un physicien, intégrer signifie la plupart du temps calculer l'aire sous la courbe de la fonction entre a et b.
La première méthode qui vienne à l'esprit, c'est de découper l'aire entre la courbe f(x), l'axe des x et les droites x= a et x = b, en une multitude de petits rectangles. Découpons l'intervalle [a,b] en rectangles élémentaires de largeur h, h étant petit. Le rectangle n° i aura donc pour longueur f(a + i*h). Sa surface est donc égale à h*f(a + i*h). L'aire sous la courbe entre a et b est obtenue en sommant tous ces petits rectangles. Reste qu'en posant cette relation, j'ai fait l'hypothèse implicite que la courbe limite le coté gauche de mon rectangle. On peut imaginer d'autres découpages. Voyons cela sur un schéma:

intégration

Comme vous le constatez, on a le choix entre trois techniques:

  • 1 - on fait coïncider le sommet haut gauche du rectangle avec la courbe : c'est la méthode des rectangles à gauche,
  • 2 - on fait coïncider le sommet haut droit du rectangle avec la courbe : c'est la méthode des rectangles à droite,
  • 3 - on fait coïncider le milieu du coté haut du rectangle avec la courbe: c'est la méthode du point milieu

Posons h = (b - a)/n, où n est le nombre de rectangles avec lesquels nous allons paver l'aire à calculer. Evidement, plus n sera grand et plus la précision du calcul sera grande (du moins en première approche!). Un rapide calcul nous montre que dans le cas:
1 - méthode des rectangles à gauche, on obtient:

\begin{align} \displaystyle \int_{a}^{b}f(x)dx \approx \sum_{i=0}^{n-1}hf(a + ih) \end{align}

2 - méthode des rectangles à droite, on obtient:

\begin{align} \displaystyle  \int_{a}^{b}f(x)dx \approx \sum_{i=0}^{n-1}hf(a + ih + \dfrac{h}{2}) \end{align}

3 - méthode du point milieu, on obtient:

\begin{align}  \displaystyle   \int_{a}^{b}f(x)dx \approx \dfrac{h}{2}(f(a) +f(b)) + h\sum_{i=1}^{n-1}hf(a + ih) \end{align}

Voilà, vous savez tout de la méthode des rectangles : très simple, facile à coder mais pas très précise. Pour des fonctions cool (polynomiales, sin, cos, exp), cette méthode donne des résultats acceptables. Et puis, on peut la programmer facilement sur une calculette...

Implémentation de la méthode des rectangles

Pour l'exemple, j'ai choisi d'implémenter la méthode du point milieu, qui est celle qui est la plus précise. j'ai donc repris la formule indiquée ci-dessu en la triturant un peu.
L'implémentation en FORTRAN (ou en C ou Python) est pratiquement immédiate. Il s'agit d'écrire une routine qui soit suffisamment générale pour être réutilisable dans tous nos programmes de physique numérique. Cela signifie qu'elle ne doit pas contenir de données propres aux programmes et que ses paramètres doivent être suffisamment complets pour supporter tous les échanges nécessaires entre le programme et la routine. Ce qui donne :

C Integration par la methode des rectangles (point milieu)
C Dominique Lefebvre janvier 2007
C
C a = borne inferieure d'integration
C b = borne superieure d'integration
C n = nombre de pas (rectangles)
C aire = surface retournee
SUBROUTINE IntRectangles (fn,a,b,n,aire)
  REAL a,b,fn,aire
  INTEGER n
  EXTERNAL fn
  REAL x,h
  C Initialisation des variables
  aire = 0
  x = a
  h = (b-a)/n
  C Boucle de calcul
  DO WHILE (x .LT. b)
    aire = aire + h*(fn(x+h)+fn(x))/2
    x = x+h
  ENDDO
END

La méthode des trapèzes

L'algorithme

La méthode des trapèzes est du même tonneau que celle des rectangles. Vous avez sans doute compris qu'on utilise non plus des rectangles pour paver l'aire mais des trapèzes. Ainsi, la partie du pavé qui jouxte la courbe est plus proche, si j'ose m'exprimer de manière aussi peu mathématicienne que cela!
Comme plus haut, je partage l'intervalle [a,b] en n petits trapèzes de largeur h = (b-a)/n. Je sais que l'aire de chaque petit trapèze est Ai = (h/2)*(f(a+ih) + f(a+(i-1)h)).
Nous obtenons l'aire recherchée en sommant l'aire de tous les trapèzes entre a et b, ce qui nous donne :

\begin{align} \displaystyle  \int_{a}^{b}f(x)dx \approx \dfrac{h}{2}(f(a) +f(b)) + h \displaystyle \sum_{i=1}^{n-1}hf(a + ih) \end{align}

La méthode des trapèzes standard est une méthode d'ordre 2, comme pourront le démontrer les fans de développements limités. Il y a toutefois une ruse pour la pousser à l'ordre 4 en estimant f"(x) par (f'(b)-f'(a))/(b-a). On appelle cette méthode la méthode des trapèzes avec correction aux extrémités. Elle donne :

\begin{align} \displaystyle  \int_{a}^{b}f(x)dx \approx \dfrac{h}{2}(f(a) +f(b)) + h\sum_{i=1}^{n-1}hf(a + ih) - \dfrac{h^2}{12}[f'(b)-f'(a)] \end{align}

Implémentation de la méthode des trapèzes en FORTRAN

L'implémentation en FORTRAN de la méthode des trapèzes standard décrit ci-dessus est pratiquement immédiate. Cela donne :

C Integration par la methode des trapezes
C Dominique Lefebvre janvier 2007
C
C a = borne inferieure d'integration
C b = borne superieure d'integration
C n = nombre de pas
C aire = surface retournee
SUBROUTINE IntTrapezes (fn,a,b,n,aire)
  REAL a,b,fn,aire
  INTEGER n
  EXTERNAL fn
  REAL h,appv INTEGER i
  C Boucle d'integration
  C Initialisation des variables
  aire = 0
  h = (b-a)/n
  C Calcul de l'aire approximative du trapeze f(a) f(b)
  app = (h/2)*(fn(a)+fn(b))
  C Boucle de calcul
  DO i=1, n-1
    aire = aire + fn(a+i*h)
  ENDDO
  C Calcul final de l'aire
  aire = app + aire*h
END

Implémentation de la méthode des trapèzes en Python

Son implémentation en Python est tout aussi immédiate :

def Trapeze(f,a,b,N):
    h   = (b-a)/float(N)
    xi  = linspace(a,b,N+1)
    fi  = f(xi)
    aire = 0.0
    for i in range(1,N):
        aire += fi[i]
    aire = (h/2)*(fi[0] + fi[N]) + h*aire
    return aire
 

Le schéma de Simpson

Dans la méthode des trapèzes, nous avons en fait interpolé f(x) par une droite entre les points i et i+h de l'intervalle. Dans la méthode de Simpson, nous n'allons plus interpoler par une droite mais par un polynôme de degré 2, ce qui va améliorer notre précision.
Plaçons nous autour d'un point x0 appartenant à l'intervalle [a,b], dans la maille de calcul x0-h et x0+h. Pour un accroissement (x-x0), le développement de Taylor limité au second ordre nous donne :
f(x) = f(x0) + (x-x0)f'(x0) + (1/2)(x-x0)^2f"(x0) + O((x-x0)^3).
Nous savons que f'(x0) = f(x0+h) - f(x0-h)/2h et que f"(x0) = (f(x0+h) - 2f(x0)+ f(x0-h))/h^2 . Si nous remplaçons ces valeurs dans le développement limité et que l'on intègre entre x0-h et x0+h, on obtient l'aire élémentaire (f(x0+h) + 4f(x0)+ f(x0-h))*h/3.
L'intégrale recherchée s'obtient en sommant toutes les aires élémentaires. Il faut quelques petites manip calculatoires sans intérêt, on obtient :

\begin{align} \displaystyle \int_{a}^{b}f(x)dx \approx \dfrac{h}{3}[(f(a) +f(b)) + 4\sum_{i=1}^{n-1}f(a + ih) + 2\sum_{i=2}^{n-2}f(a + ih)] \end{align}

La méthode de Simpson est une méthode d'ordre 4. Vous remarquerez l'alternance des sommes en facteur de 2 et de 4. On se sert de cette caractéristique pour simplifier la programmation.
ATTENTION : le nombre d'intervalles doit être pair.

Implémentation du schéma de Simpson en C

Le schéma de Simpson peut être implémenté de la façon suivante:

C Integration par la methode de Simpson
C Dominique Lefebvre janvier 2007
C
C a = borne inferieure d'integration
C b = borne superieure d'integration
C n = nombre de pas
C aire = surface retournee
SUBROUTINE IntSimpson(fn,a,b,n,aire)
  REAL a,b,fn,aire
  INTEGER n
  EXTERNAL fn
  REAL h,SommePaire, SommeImpaire
  INTEGER i
  C Boucle d'integration
  C Initialisation des variables
  aire = 0
  h = (b-a)/(n*2)
  SommePaire = 0
  SommeImpaire = 0
  C Calcul de la somme des indices impaires
  DO i=1, n-1
   SommeImpaire = SommeImpaire + fn(a+h*2*i)
  ENDDO
  C Calcul de la somme des indices paires
  DO i=1, n
    SommePaire = SommePaire + fn(a+h*(2*i-1))
  ENDDO
  C Calcul final de l'aire
  aire = h*(fn(a) + fn(b)+ 2*SommePaire + 4*SommeImpaire)/3
END

Implémentation du schéma de Simpson en Python

Voici l'implémentation en Python de la fonction Simpson() :

def Simpson(f,a,b,n):
    if n % 2 != 0 :
        return -1
    else :       
        h=(b-a)/float(n)
        sum1=(f(a)+f(b))/6.0
        for i in range(1,n) :
           sum1 += f(a+i*h)/3.0
        for i in range(n) :
           sum1 += f(a+(2*i+1)*h/2.0)*2.0/3.0
    I = h*sum1
    return I

Vous noterz que j'ai ajouté ici un contrôle de la parité du nombre d'intervalles. Si ce nombre n'est pas pair, la fonction retourne la valeur -1. Vous trouverez un exemple d'utilisation de cette fonction dans le script ExempleSimpson.py.

Le schéma de Romberg

L'algorithme

Le schéma de Romberg utilise une extrapolation de Richardson à partir de la méthode des trapèzes. Pour plus d'explications, je vous renvoie à votre cours d'analyse numérique ou à l'article de Wikipedia Méthode_de_Romberg dont je me suis inspiré.

Implémentation du schéma de Romberg en C

Pour coder ce schéma en C, j'ai emprunté le principe de calcul en tableau à l'excellent bouquin de A. Garcia, que j'ai maintes fois cité : "Numerical Methods for Physics". La programmation est un peu moins évidente que celle des méthodes précédentes. Il y a plus élégant, aussi ! Par exemple, la programmation de cette méthode dans les "Numerical Recipes in Fortran" est plus fine mais bon ! La routine ci-dessous fonctionne aussi...

C Integration par la methode de Romberg
C Dominique Lefebvre fevrier 2007
C Inspire du code de A.Garcia dans Numerical Methods for Physics
C
C Modifiee le 30/08/09 suite à une erreur d'estimation de la precision
C signalee par Melle Nodet (Universite de Grenoble)
C
C a = borne inferieure d'integration
C b = borne superieure d'integration
C eps = précision souhaitee
C aire = surface retournee
SUBROUTINE IntRomberg(fn,a,b,eps,aire)
  REAL fn,a,b,eps,aire
  EXTERNAL fn
  INTEGER out, pieces, i, nx(16)
  INTEGER nt, ii, n, nn, l, ntra, k, m, j, l2
  REAL h,xi, delta, sum
  REAL c, fotom, t(136)
  LOGICAL fini
  C Calcul du premier terme
  h = b-a
  pieces = 1
  nx(1) = 1
  delta = h / pieces
  c = (fn(a) + fn(b)) / 2.0
  sum = c
  t(1) = delta * c
  n = 1
  nn = 2
  fini =.false.
  DO WHILE (.not. fini)
    n = n + 1
    fotom = 4.0
    nx(n) = nn
    pieces = pieces * 2
    l = pieces - 1
    l2 = (l+1) / 2
    delta = h / pieces
    C Evaluation par la methode des trapezes
    DO ii = 1, l2
      i = ii * 2 - 1
      xi = a + delta * i
      sum = sum + fn(xi)
    ENDDO
    t(nn) = sum * delta
    ntra = nx(n-1)
    k = n-1
    C Calcul de la nieme colonne du tableau
    DO m = 1, k
      j = nn + m
      nt = nx(n-1) + m - 1
      t(j) = (fotom * t(j-1) - t(nt)) / (fotom - 1)
      fotom = fotom * 4
    ENDDO
    C Estimation de la precision
    IF (n .lt. 5) THEN
      nn = j+1
    ELSE IF (t(nn+1) .eq. 0.0) THEN
      nn = j+1
    ELSE
     IF
      (abs(t(ntra+1)-t(nn+1)) .le. abs(t(nn+1)*eps)) THEN
      fini = .true.
     ELSEIF (abs(t(nn-1)
      - t(j)) .le. abs(t(j)*eps)) THEN
      fini = .true.
     ELSE
      nn = j+1
     ENDIF
    ENDIF
  ENDDO
  C On retourne la valeur finale de l'aire
  aire = t(j)
  RETURN
END

Implémentation du schéma de Romberg en Python

En Python, j'ai utilisé le même algorithme mais en l'adaptant aux possibilités de Python. Cela nous donne un code beaucoup plus compact. J'utilise la fonction Trapeze() décrite plus haut :

def Romberg(fonction,a,b,eps,nmax):
    A = zeros((nmax,nmax),float)
    for i in range(0,nmax):
        N = 2**i
        A[i,0] = Trapeze(fonction,a,b,N)
        for k in range(0,i):
            n = k + 2
            A[i,k+1] = 1.0/(4**(n-1)-1)*(4**(n-1)*A[i,k] - A[i-1,k])
        if (i > 0):
            if (abs(A[i,k+1] - A[i,k]) < eps):
               break  
    return A[i,k+1]

Vous trouverez un exemple d'utilisation de cette fonction dans le script ExempleRomberg.py.

Un programme pour tester les différentes méthodes d'intégration

J'ai fait un petit programme FORTRAN qui permet, sur une fonction simple dont on connait le résultat, de comparer la précision des quatre méthodes. Je définis d'abord la fonction à intégrer, ici x^2, dont la primitive est bien connue. La fonction à intégrer est définie dans la routine fn(x) décrite ci-dessous:

C Definition de la fonction a integrer
C Dominique Lefebvre Janvier 2007
REAL FUNCTION fn(x)
  REAL x
  fn = x*x
  RETURN
END

Le corps du programme d'intégration, qui appelle les quatre schémas d'intégration est le suivant:

PROGRAM TestIntegration
C Programme de test des differentes methodes d'integration
C des fonctions reelles.
C
C Dominique Lefebvre janvier 2007
IMPLICIT NONE
C Declaration de la fonction a integrer
REAL fn
EXTERNAL fn
C Declaration des subroutines d'integration
EXTERNAL IntRectangles
EXTERNAL IntTrapezes
EXTERNAL IntSimpson
EXTERNAL IntRomberg
C Declaration des variables
REAL a,b,s,s1,eps
INTEGER hinit
C Saisie des parametresv WRITE(*,'(a,$)') 'Borne inferieure d''integration: '
READ (*,*) a
WRITE(*,'(a,$)') 'Borne superieure d''integration: '
READ (*,*) b
WRITE(*,'(a,$)') 'Precision requise: '
READ (*,*) eps
C Premier calcul avec un pas = 4096v hinit = 4096
CALL IntRectangles(fn,a,b,hinit,s1)
C Initialisation de s
s = s1 + 2*eps
C Boucle pour atteindre la precision voulue
DO WHILE (ABS(s-s1) .GE. eps)
hinit = hinit*2
CALL IntRectangles(fn,a,b,hinit,s)
ENDDO
WRITE(*,*)'La valeur de l''integrale (Rectangles) est : ', s
C Calcul par la methode des trapezes
hinit = 4096
s = 0
CALL IntTrapezes(fn,a,b,hinit,s)
WRITE(*,*)'La valeur de l''integrale (Trapezes) est : ', s
C Calcul par la methode de Simpson
hinit = 8096
s = 0
CALL IntSimpson(fn,a,b,hinit,s)
WRITE(*,*)'La valeur de l''integrale (Simpson) est : ', s
C Calcul par la methode de Romberg
s = 0
CALL IntRomberg(fn,a,b,eps,s)
WRITE(*,*)'La valeur de l''integrale (Romberg) est : ', s
STOP
END

Vous pouvez maintenant exécuter le programme, par exemple pour intégrer x2 entre 0 et 1. Le résultat est connu analytiquement, c'est 1/3. Voyons ce que donne l'exécution :

I:\NumLab\LabPhysique\Integration\Output>integration
Borne inferieure d'integration: 0
Borne superieure d'integration: 1
Precision requise: 1.0e-5
La valeur de l'integrale (Rectangles) est : 0.333333969
La valeur de l'integrale (Trapezes) est : 0.333333552
La valeur de l'integrale (Simpson) est : 0.333312601
La valeur de l'integrale (Romberg) est : 0.333333343

On constate sans surprise que la méthode de Romberg donne les meilleurs résultats. Mais les autres méthodes, en particulier les trapèzes, sont tout à fait convenables pour les besoins courants.

Tous les codes de cette page sont téléchargeables dans la bibliothèque de codes de TangenteX.

Contenu et design par Dominique Lefebvre - www.tangenteX.com janvier 2007   Licence Creative Commons   Contact :