Méthode de Romberg

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

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, du type $[a+kh,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+kh)\right]$

où les pondérations $\omega _{k}$ sont égales à 1 pour les points extrêmes, et à 2 pour les autres. La méthode est dite d’ordre 1, car donne un résultat exact pour un polynôme de degré 1. D'après la formule d'Euler-Maclaurin, 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 .$

où les coefficients c_k ne dépendent pas du pas h choisi, et seules les puissances paires de h apparaissent. Cette propriété (qui est aussi partagée avec la méthode du point médian) sera avantageusement exploitée, pour supprimer, un à un les termes de l'erreur, à l'aide de l'extrapolation de Richardson.

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

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^{m}R(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

Critère d’arrêt

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

Redondances

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ⁿ subdivisions pour des valeurs croissantes de n.

D'après la formule suivante :

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

avec $hk=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

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

On considère le calcul de l'intégrale

$I=\int _{-1}^{2}{\frac {2}{1+4t^{2}}}\,\mathrm {d} t.$

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

• Méthode des trapèzes
• Méthode de Simpson
• Extrapolation de Richardson

Lien externe

• (en) Recursive integration formulas from Romberg integration

.