Méthode de Romberg

Un article de Wikipédia, l'encyclopédie libre.
Aller à : navigation, rechercher

En analyse numérique, la méthode d'intégration de Romberg est une méthode récursive de calcul numérique d'intégrale, fondée sur l'application du procédé d'extrapolation de Richardson à la méthode des trapèzes. Cette technique d'accélération permet d'améliorer l'ordre de convergence de la méthode des trapèzes, en appliquant cette dernière à des divisions dyadiques successives de l'intervalle d'étude et en formant une combinaison judicieuse.

Principe[modifier | modifier le code]

On souhaite calculer l'intégrale I = \int_a^b f(x) \, dx d'une fonction f supposée régulière sur [a,b] en subdivisant [a,b] en n sous-intervalles identiques (n pair, n = 2^p par exemple), du type [a + k h, a + (k+1) h] pour k = 0, ..., n-1 et h=(b-a)/n. La méthode des trapèzes, notée T(n), est définie à l'aide de cette grille régulière :

T(n) = \frac{h}{2} \left[\sum_{k=0}^{n} \omega_k f(a + k h)\right]

où les pondérations \omega_k sont égales à 1 pour les points extrêmes, et à 2 pour les autres. Puisque la méthode est exacte pour un polynôme de degré 1 (méthode d’ordre 1), l'erreur commise, notée E = I-T(n)\,, vérifie

E(h) =  c_1 h^2 + c_2 h^4 + c_3 h^6 \cdots.

Cette relation exprime le fait que la méthode des trapèzes présente une erreur proportionnelle à h^2 ; on dit qu'elle est en O(h^2).

Des manipulations algébriques fournissent une approximation de I en O(h^4). Pour ce faire, il s'agit d'éliminer le terme en h^2 dans le développement de l'erreur. On y parvient en considérant la quantité [4 T(n) – T(n/2)]/3, qui correspond précisément à la méthode de Simpson. Le même procédé peut être reconduit, permettant ainsi d'annuler successivement les termes en h^4, h^8, \ldots, conduisant ainsi à des approximations de I qui sont de plus en plus précises.

Algorithme[modifier | modifier le code]

Formalisations de la technique précédente de réduction de l'erreur :

  • Initialisations :
R(0,0) = \frac{1}{2} (b-a) (f(a) + f(b)).
R(n,0) = T(2^n) , \,
soit le résultat de la méthode des trapèzes basée sur
 h_n = \frac{b-a}{2^n}.
  • Récurrence :
R(n,m) = \frac{4^mR(n, m-1)-R(n-1,m-1)}{4^{m}-1}

On montre par récurrence que l'approximation à l’étape n, soit R(n,n), fournit une approximation de I qui est en O(h_n^{2n+2}), ceci à condition que la fonction soit de classe \mathcal{C}^{2n+2} (ou \mathcal{C}_I^{2n+2}). Le calcul de R(n,n) est exact pour les polynômes de degré 2 n + 1 : elle est d’ordre 2n + 1.

L'algorithme peut être représenté sous la forme d'un tableau. La première diagonale R(n,n) fournit les approximations successives ; la première colonne R(n, 0) correspond (par définition) à la méthode du trapèze, la seconde R(1,n) reflète la méthode de Simpson, la troisième celle de Boole qui est la formule de Newton-Cotes à 5 points appliquée à une décomposition en 2^{n-2}+1 sous-intervalles. Par contre, les colonnes suivantes correspondent à des formules de quadrature différentes de celles de Newton-Cotes d’ordre supérieur, évitant ainsi les problèmes d’instabilité.

En pratique[modifier | modifier le code]

Critère d’arrêt[modifier | modifier le code]

Pour obtenir une approximation de I avec une précision \epsilon>0 donnée, le critère d'arrêt est satisfait lorsque |R(n, n)-R(n-1, n-1)|<\epsilon. Ceci implique généralement |R(n, n)-I|<\epsilon.

Redondances[modifier | modifier le code]

Lorsque l’évaluation de la fonction f (en un point) comporte un certain coût de traitement, il est préférable d’éviter d’appliquer brutalement la méthode du trapèze avec 2^n\, subdivisions pour des valeurs croissantes de n.

D'après la formule suivante :

2 T(2 k) = T(k) + h \sum_{j=1}^k f(a + (j-1/2)h)

avec h k = b - a, l’initialisation R(n,0) = T(2^n)\, est avantageusement remplacée par la relation récursive :

R(n,0) = \frac{1}{2} R(n-1,0) + h_n \sum_{j=1}^{2^{n-1}} f(a + (2j-1)h_n).

Ceci permet de réduire d’un facteur 2 le nombre total d’évaluations de f qui, pour déterminer R(n,0), atteint ainsi 2^n+1.

Subdivision initiale[modifier | modifier le code]

Décomposer [a, b] en p morceaux avant d’appliquer la méthode de Romberg sur chacun d’eux peut être utile lorsque la fonction n’est pas assez régulière (limitation de n), ou encore lorsqu’elle est régulière uniquement sur chacun des morceaux.

Dans le cas contraire, une subdivision initiale présente généralement peu d’intérêt. Il ne faut cependant pas s'attendre à une amélioration considérable de la précision par la seule augmentation du nombre d'itérations : en effet, bien que le terme h_n^{2n+2} diminue à chaque itération, le coefficient multiplicateur peut augmenter considérablement puisqu'il se réfère à des dérivées d'ordres plus élevés.

Exemple[modifier | modifier le code]

A intégrer f(t) = 2/(1 + 4t^2)\, sur [-1,2]. On trouve I = \arctan(4)+\arctan(2) \simeq 2,432966381462123. La matrice de l'algorithme est


k=0 1 2 3 4
p=0 0,7764705882



1 1,8882352941 2,2588235294


2 2,3510141988 2,5052738337 2,5217038540

3 2,4235526286 2,4477321053 2,4438959900 2,4426609446
4 2,4307735880 2,4331805744 2,4322104723 2,4320249879 2,4319832783

d'où l'on tire les approximations consécutives

k R(k)
0 0,77647058823529
1 2,25882352941176
2 2,52170385395538
3 2,44266094457555
4 2,43198327829659

Voir aussi[modifier | modifier le code]

Lien externe[modifier | modifier le code]