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:

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:

intégration

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

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:

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.


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

Cette œuvre est mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Pas de Modification 3.0 France.