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 voire simplistes. 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, on en reparlera...

La méthode des rectangles

Considérons donc une fonction continue sur un intervalle [a,b]. Je ne vais pas vous faire un cours sur l'intégration! Pour un physicien, intégrer signifie 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 de largeur faible, appelons la h, et de hauteur f(h).

L'aire sous la courbe est obtenue en sommant tous ces petits rectangles. Voyons cela sur un schéma:

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

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

3 - méthode du point milieux, on obtient

Voilà, vous savez tout de la méthode des rectangles: très simple mais pas très précise. Mais facile à coder! Pour des fonctions cool (polynomiales, sin, cos, exp), cette méthode donne des résulats 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 est pratiquement immédiate. Il s'agit d'écrire une routine qui soit suffisamment générale pour être réutilisable dans tous nos programme 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

Le programme colle à l'algorithme. Il n'y a pas de piège caché...

HOME

La méthode des trapèzes

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 :

La méthode des trapèzes standard est une méthode d'ordre 2, comme pourront le démontrer les fans du développement de Taylor. 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 :

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

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

Ce qui 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,app
    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

Là non plus, ce code ne nécessite pas vraiment de commentaire.

HOME

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 (vous voulez vraiment savoir? Si oui, passez moi un mail ou bien regardez dans un cours d'analyse numérique). Et l'on obtient :

La méthode de Simpson est une méthode d'ordre 4.

ATTENTION : le nombre d'intervalles doit être pair.

Vous remarquerez l'alternance des sommes en facteur de 2 et de 4. On se sert de cette caractéristique pour simplifier la programmation. 

Implémentation du schéma de Simpson

Ce qui donne:

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

HOME

Le schéma de Romberg

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 dont je me suis inspiré http://fr.wikipedia.org/wiki/Methode_de_Romberg.

Pour coder ce schéma, 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": un excellent investissement!

Implémentation du schéma de Romberg

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...