Les méthodes d'intégration numérique
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 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, voir par exemple la méthode de Monte-Carlo...
La méthode des rectangles
L'algorithme
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. 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:
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} \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} \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 milieux, on obtient:
\begin{align} \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 mais pas très précise. Mais facile à coder! 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
Le programme colle à l'algorithme. Il n'y a pas de piège caché...
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} \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}
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 :
\begin{align} \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
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,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.
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 :
\begin{align} \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.
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
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
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 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...
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 Maelle 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
Un programme pour tester les différentes méthodes d'intégration
J'ai fait un petit programme 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 (oui x^3/3 + constante..). 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 parametres
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 = 4096
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
Comme d'habitude, il faut saisir le programme dans VFort, le compiler et l'exécuter. je ne reviens plus sur les manips...
Pour les paresseux, voici les codes:
- La fonction à intégrer: IntFonction.for
- Les différents schémas: IntRectangles.for , IntTrapezes.for , IntSimpson.for, IntRomberg.for
- Le programme de test: TestIntegration.for
- Le projet VFort : Integration.prj
Vous pouvez maintenant exécuter le programme, par exemple pour intégrer x^2 entre 0 et 1. Le résultat est connu analytiquement, c'est 1/3 bien sur!
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.
Contact : PhysiqueX ou 