TangenteX.com
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 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...
Considérons donc une fonction de \( \mathbb{R} \) dans \( \mathbb{R}
\) continue sur un intervalle [a,b]. 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:
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...
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 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}
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
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
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.
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
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 ExempleSimpson.py.
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 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
En Python, j'ai utilisé le même algorithme mais en l'adaptant aux
possibilités de Python. Cela nous donne un code 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]
Quande je disais plus compact, je n'exagérais pas ! Vous trouverez un exemple d'utilisation de cette fonction dans le script ExempleRomberg.py.
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
Comme d'habitude, il faut saisir le programme dans VFort, le compiler et l'exécuter. je ne reviens plus sur les manips. Voici les codes téléchargeables :
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.
Vous pourrez construire le même programme en Python à partir des exemples que j'ai cité plus haut.