- Calcul numérique d'une intégrale
-
En analyse numérique, il existe une vaste famille d’algorithmes dont le but principal est d’estimer la valeur numérique de l’intégrale définie sur un domaine particulier pour une fonction donnée (par exemple l’intégrale d’une fonction d’une variable sur un intervalle).
Ces techniques procèdent en trois phases distinctes :
- Décomposition du domaine en morceaux (un intervalle en sous-intervalles contigus).
- Intégration approchée de la fonction sur chaque morceau.
- Sommation des résultats numériques ainsi obtenus.
On appelle formule de quadrature une expression linéaire dont l’évaluation fournit une valeur approchée de l’intégrale sur un morceau typique (l’intervalle [0, 1] par exemple). Une transformation affine permet de transposer la formule sur un morceau particulier. La formule de quadrature fait intervenir des valeurs pondérées de la fonction (et parfois également celles de sa dérivée) en certains nœuds : les coefficients de pondération et les nœuds dépendent de la méthode employée. Ces formules de quadrature sont en effet obtenues à l’aide de la substitution de la fonction par une approximation, c’est-à-dire par une fonction proche dont l’intégrale peut être déterminée algébriquement.
Une indication grossière de l’efficacité d’une formule de quadrature est son ordre qui, par définition, est la plus grande valeur entière m pour laquelle la valeur approchée de l’intégrale soit exacte pour tout polynôme de degré inférieur ou égal à m.Cependant, la précision du résultat obtenu dépend à la fois de l’ordre de la formule de quadrature, de la taille des morceaux et de la régularité de la fonction. D’autre part, il est généralement inutile d’appliquer une formule de quadrature d’ordre m si la fonction n’est continûment dérivable jusqu’à l’ordre m + 1.
Sommaire
- 1 Méthode de calcul d'intégrale à une dimension
- 2 Autres méthodes de quadrature numérique
- 3 Notes (et références)
- 4 Voir aussi
Méthode de calcul d'intégrale à une dimension
Généralités
Considérons une intégrale définie dont on cherche à estimer la valeur numérique.
Hypothèses et traitements préalables
Supposons que a et b soient finis : dans le cas contraire, il est conseillé d’effectuer un changement de variable permettant de satisfaire cette hypothèse[1].
Supposons également que la fonction f à intégrer ne comporte pas de singularité. Par exemple, la fonction f(x) = x − α avec 0 < α < 1 est intégrable sur [0,1] et I = 1 / (1 − α). Bien que f soit parfaitement régulière sur (0,1], la singularité en 0 et l’impossibilité de la prolonger par une fonction continue cause de grandes difficultés à toutes les méthodes d’intégration numérique, en particulier pour celles qui utilisent explicitement f(0)[2].Dans une telle situation, il convient de soustraire à f une fonction g dont l’intégrale est connue et qui soit telle que f − g ne soit plus singulière, puis d’intégrer numériquement cette différence.
Mise en œuvre
Considérons une formule de quadrature associée à [0,1] du type où les pondérations αi et les nœuds xi sont donnés.
Partant d’une décomposition régulière de en n sous- intervalles de longueur h = (b − a) / n, soit les intervalles pour l’application de la formule de quadrature précédente à chaque Jk s’effectue à l’aide d’une transformation affine, permettant ainsi d’obtenir une approximation In(f) de I qui s’écrit :
Cette relation est la formule composite associée à une formule de quadrature générale. L’ordre du calcul des termes de cette double somme et certains arrangements permettent le plus souvent de réduire le nombre d’opérations (évaluations de f(.) et multiplications par αi). Cette question est développée plus loin pour quelques formules de quadrature particulières.
Si la formule de quadrature comporte des termes du type faisant intervenir la dérivée ou du type impliquant la dérivée seconde, la transformation affine fait apparaître des facteurs h ou h2, ceci conformément à la relation suivante :
Ordre de la formule de quadrature et convergence
Si m est l’ordre de la formule de quadrature et si f(x) est de classe (soit l’espace des fonctions m + 1 fois dérivables dont la dérivée m + 1 est continue par morceaux), notons
- .
Dans ce cas, il existe une constante C indépendante de f et de telle que
PreuvePour tout la substitution de f(x) restreinte à Jk par son développement de Taylor d’ordre m autour de la borne inférieure x0,k = a + kh implique que le remplacement de f(x) par le reste Rm,k(x) n'affecte pas la valeur de l’écart : l’ordre m assure en effet l’exactitude de la formule de quadrature pour chaque terme polynomial.
Par ailleurs, l’inégalité de Taylor-Lagrange appliquée à chaque Jk implique :
- pour tout
Il en découle
En choisissant il vient
Remarque : ce développement n’a pas pour objectif de déterminer la constante C la plus faible.Ce résultat conforte les recommandations suivantes :
- Si fm(x) n’est pas continue sur , une formule de quadrature d’ordre m (ou plus) présente peu d’intérêt.
- Si f(x) est régulière par morceaux, il vaut la peine de décomposer en sous-intervalles correspondant aux morceaux de régularité, puis d’appliquer une formule composite de quadrature sur chaque morceau.
- La même approche peut se révéler opportune pour intégrer une fonction régulière, mais dont la variabilité est très dissemblable d’une zone à l’autre. L’intérêt se manifeste toutefois principalement sur le volume des calculs.
Formules simples
Ces méthodes utilisent l’interpolation des fonctions à intégrer par des polynômes dont la primitive est connue.
Formules du rectangle et du point milieu
C’est la méthode la plus simple qui consiste à interpoler la fonction f à intégrer par une fonction constante (polynôme de degré 0).
Si ξ est le point d’interpolation, la formule est la suivante :
Le choix de ξ influence l’erreur :
- Si ou , l’erreur est donnée par
- .
- C’est la méthode du rectangle qui est d’ordre 0.
- Si , l’erreur est donnée par
- .
- Il s’agit de la méthode du point milieu qui est d’ordre 1.
Ainsi, le choix du point milieu améliore l’ordre de la méthode : celle du rectangle est exacte (c’est-à-dire E(f) = 0) pour les fonctions constantes alors que celle du point milieu est exacte pour les polynômes de degré 1. Ceci s’explique par le fait que l’écart d’intégration de la méthode du point milieu donne lieu à deux erreurs d’évaluation, de valeurs absolues environ égales et de signes opposés.Formule du trapèze
En interpolant f par un polynôme de degré 1, les deux points d'interpolation et suffisent à tracer un segment dont l’intégrale correspond à l’aire d’un trapèze, justifiant le nom de méthode des trapèzes qui est d’ordre 1 :
L’erreur est donnée par
Conformément aux expressions de l’erreur, la méthode des trapèzes est souvent moins performante que celle du point milieu.
Formule de Simpson
En interpolant f par un polynôme de degré 2 (3 degrés de liberté), 3 points (ou conditions) sont nécessaires pour le caractériser : les valeurs aux extrémités a, b, et celle choisie en leur milieu x1 / 2 = (a + b) / 2. La méthode de Simpson est basée sur un polynôme de degré 2 (intégrale d’une parabole), tout en restant exacte pour des polynômes de degré 3 ; elle est donc d’ordre 3 :
L’erreur globale est donnée par
Remarque : comme la méthode du point milieu qui caractérise un polynôme de degré 0 et qui reste exacte pour tout polynôme de degré 1, la méthode de Simpson caractérise un polynôme de degré 2 et reste exacte pour tout polynôme de degré 3. Il s’agit d’une sorte d’anomalie où se produisent des compensations bénéfiques à l’ordre de la méthode.
Généralisation
La fonction f peut être interpolée à l’aide de son évaluation en m points équidistants (comprenant les deux extrémités si m > 1, méthode du point milieu si m = 1) par un polynôme de degré m − 1 issu d’une base de polynômes de Lagrange et dont l’intégrale est fournie par les formules de Newton-Cotes. Ce procédé permet ainsi une généralisation des résultats précédents. Avec m points, il en découle une méthode
- d’ordre m si m est impair (anomalie),
- d’ordre m − 1 si m est pair.
On notera ici NC-m la méthode basée sur m points :
- NC-1 est la méthode du point milieu
- NC-2 est la formule du trapèze
- NC-3 est la formule de Simpson
Pour des questions d’instabilité numérique provenant en particulier du phénomène de Runge, il est cependant préférable de limiter le degré m du polynôme d'interpolation, quitte à subdiviser l'intervalle en sous-intervalles.
Formule faisant intervenir la dérivée
Formule NC-m-m'
Il s’agit d’une généralisation des formules NC-m dans lesquelles interviennent non seulement la fonction évaluée en m points équidistants, mais également la dérivée de la fonction évaluée en m' points équidistants ; malgré l’abus de langage, on notera ici NC-m-m' une telle formule.
On se limitera ici à m' = 2 correspondant aux deux extrémités a et b.
Peu connues (et donc rarement présentées dans les cours), ces méthodes permettent de gagner deux ordres de convergence par rapport à la méthode correspondante sans la dérivée, ceci en nécessitant très peu de calculs supplémentaires : en effet, les coefficients de f'(a) et de f'(b) sont opposés et ainsi, dans la formule composite (dont il est question ci-dessous), les dérivées aux extrémités de deux intervalles adjacent se simplifient.
Si xi désigne les points d’évaluation de f (i entre 0 et m − 1) :
- Formule NC-1-2 : basée sur un polynôme de degré 2, elle est d’ordre 3 :
- Formule NC-2-2 : basée sur un polynôme de degré 3, elle est d’ordre 3 :
- Formule NC-3-2 : basée sur un polynôme de degré 4, elle est d’ordre 5 :
- Formule NC-4-2 : basée sur un polynôme de degré 5, elle est d’ordre 5 :
- Formule NC-5-2 : basée sur un polynôme de degré 6, elle est d’ordre 7 :
Concernant l’erreur globale d’une formule de quadrature linéaire d’ordre p, elle est donnée par
Ce procédé permet ainsi une généralisation des formules de Newton-Cotes. Avec m points pour la fonction et 2 points pour sa dérivée, il en découple une méthode- d’ordre m + 2 si m est impair (anomalie),
- d’ordre m + 1 si m est pair.
Formule NC-m-m'-m"
La notation indique que la dérivée seconde intervient également en m'' points équidistants. Mentionnons uniquement une formule particulièrement remarquable qui présente une double anomalie :
- Formules NC-4-2-2 : basée sur un polynôme de degré 7, elle est d’ordre 9 :
- Les coefficients de f'(a) et de f'(b) sont opposés, ce qui permet des simplifications dans la formule composite. Par contre, ce n’est pas le cas pour f''.
Remarques :
- Une méthode NC-1-1-...-1 n’est autre que l’intégration d’un développement de Taylor.
- Dans une formule NC-m-m'-m", m doit être positif afin de tenir compte d’une translation sur f'.
- Tout triplet d’entiers (m, m', m") ne conduit pas nécessairement à une formule de quadrature. C’est le cas pour les triplets du type (2 p, 2 q +1, r), ainsi que quelques cas du type (2 p + 1, 2 q, 2 r +1) [3].
- Les triplets du type (2 p +1, 2 q, 2 r), (2 p +1, 2 q + 1, 2 r + 1) et (2 p, 2 q, 2 r +1) présentent une anomalie dans le sens où la formule est basée sur un polynôme de degré m + m' + m" - 1 (nombre de degrés de liberté) alors qu’elle est encore vraie pour un polynôme de degré m + m' + m".
Formules composites
Pour chacune des méthodes précédentes, le terme d’erreur dépend d’une puissance de b-a. Cette imprécision étant le plus souvent trop importante, l’erreur peut être réduite en découpant simplement l’intervalle [a, b] en n sous-intervalles (supposés de longueurs égales), ceci dans le but de déterminer une valeur approchée de l’intégrale sur chacun d’eux, en application de la méthode choisie. L’intégrale sur [a, b] est estimée par la somme des valeurs ainsi calculées.
On appelle formule composite l’expression caractérisant cette estimation.
Notons k l’indice des n sous-intervalles, h = (b − a) / n la longueur de chacun d’eux, xk = a + kh la borne inférieure et mk = a + (k + 1 / 2)h le point milieu, ceci pour k entre 0 et n − 1. La formule composite sont les suivantes :
- Méthode du point milieu d’ordre 1 :
- Méthode des trapèzes d’ordre 1 :
- Méthode de Simpson d’ordre 3 :
Pour les formules de quadrature faisant intervenir des dérivées, donnons deux exemples, les autres cas se déduisant par analogie (ici xk,i = a + kh + iΔ où Δ = h / (m − 1), m étant le nombre de points d’un sous-intervalle) :
- NC-3-2 d’ordre 5 :
- NC-4-2-2 d’ordre 9 :
Concernant l’erreur finale d’une formule de quadrature linéaire d’ordre p, elle est donnée par la relation
Méthode de Romberg
Article détaillé : méthode de Romberg.S’agissant de choisir une méthode d’intégration numérique, cette méthode ne doit pas être négligée, en particulier lorsque la fonction est régulière. Après n itérations, elle conduit à une méthode d’ordre de 22n + 1 avec une erreur en O(h2n + 2) pour des fonctions de classe , ceci avec une grande économie du nombre d’évaluations de la fonction (précisément 2n + 1 évaluations).
Tirages aléatoires
Article détaillé : Méthode de Monte-Carlo.Pour intégrer une fonction f sur un intervalle [a,b], la méthode de Monte-Carlo est ici mentionnée à titre "presque anecdotique" : sa performance reste en effet très limitée et son coût de traitement élevé à cause du grand nombre d’évaluations de f qui sont nécessaires pour espérer obtenir un résultat significatif.
L’estimation de l’intégrale I de f est fournie par Iq(f) = (b − a)Mq(f)) où Mq(f)) est la moyenne arithmétique des f(ui), les ui étant q tirages aléatoires indépendants et distribués uniformément sur [a,b].
Cette méthode est d’ordre 0 puisqu’elle donne un résultat exact pour toute fonction constante (même avec un unique tirage).
Cependant, ne s’agissant pas d’une formule de quadrature, l’erreur ne peut pas être majorée avec certitude par une quantité décroissante avec le nombre de tirages. En fait, à chaque groupe de q tirages correspond une estimation particulière de l’intégrale, c'est-à-dire une réalisation d’une variable aléatoire Iq(f) dont la distribution dépend de f, de [a,b] et de q. On peut alors assurer que son espérance est égale à l’intégrale de f et préciser une borne de l’écart type de l’erreur.
Les résultats suivants permettent de caractériser la distribution de l’erreur et son écart type :- Lorsque q augmente, l’erreur Eq(f) = I − Iq(f) s’amenuise statistiquement. Plus précisément :
- Soit σ l’écart type d’un tirage qui est défini par la relation
- Alors la distribution de converge vers celle de la loi normale centrée réduite, soit . En d’autres termes, Eq(f) est essentiellement distribué comme une loi normale
- Si f est une fonction bornée quelconque (la borne est notée ), alors et ainsi
- .
- Si f est de classe (dérivable et à dérivée continue, en particulier bornée[4]), alors et ainsi
- Si f est de classe , que [a,b] est décomposé en n intervalles de longueur h avec nh = b − a, que le nombre total q de tirages est réparti à parts égales dans chaque intervalle, alors
PreuveLes tirages :
Soit u une variable aléatoire de loi uniforme sur [a,b] et de densité
La variable aléatoire f(u) est alors d’espérance
et de variance
1er résultat :C’est une conséquence directe du théorème central limite.
2e résultat :La conclusion découle alors du 1er résultat.
Cette inégalité peut être montrée plus rigoureusement à l’aide de l’inégalité de Hölder pour tout
3e résultat :Puisque la méthode de Monte-Carlo est d’ordre 0, on peut supposer sans restreindre la généralité que l’intégrale de f est nulle. Dans ce cas, le théorème de Rolle appliqué à une primitive de f implique l’existence d’un point tel que f(μ) = 0. A l’aide du théorème de Taylor, pour tout , . Par conséquent,
d’où le résultat.
4e résultat :Notons Ji l’un des intervalles, puis gi(u) la fonction f(u) restreinte à Ji après soustraction d’une constante égale à la moyenne de f sur Ji. Ainsi gi'(u) = f'(u) et l’intégrale de gi sur Ji est nulle.
On utilise maintenant le 3e résultat pour caractériser l’erreur d’intégration de f sur Ji issue de qi tirages, soit qui, puisque la méthode de Monte-Carlo est d’ordre 0, est égale à :
En utilisant l’indépendance des variables et l’hypothèse que les qi sont tous égaux à q / n, il vient
Remarque :
Le dernier résultat justifie le comportement relativement médiocre de la méthode de Monte-Carlo. Supposons en effet que le nombre q de tirages soit imposé (limitation de l’effort de traitement). Dans ce cas et afin de réduire σ(Eq(f)), il convient de choisir le plus petit h possible, soit un tirage par intervalle (n = q) pour obtenir
Autant choisir la méthode du point milieu si la fonction est dans .
Tableau comparatif
Le tableau suivant résume les performances théoriques de chaque méthode :
Nom de la
MéthodeDegré du
polynômeNombre de
pointsDegré
d’exactitude
(ordre)Degré
d’erreur globaleDegré
d’erreur finaleRectangle 0 1 0 2 1 Point Milieu 0 1 1 3 2 Trapèze 1 2 1 3 2 Simpson 2 3 3 5 4 NC-5 4 5 5 7 6 NC-1-2 2 1+2 3 5 4 NC-3-2 4 3+2 5 7 6 NC-5-2 6 5+2 7 9 8 NC-4-2-2 7 4+2+2 9 11 10 Romberg (n) 2n+1 2n+1 2n+1 - 2n+2
Signification des titres des colonnes :- Degré du polynôme : degré des polynômes sur lesquels se base la formule.
- Nombre de points : sur chacun desquels la fonction est évaluée.
- Degré d’exactitude : degré maximal des polynômes pour lesquels la formule est exacte (c’est l’ordre de la méthode).
- Degré d’erreur globale : la puissance du facteur b − a dans l’erreur pour une application globale de la formule sur l’intervalle.
- Degré d’erreur finale : la puissance du facteur h dans la formule composite.
Exemple numérique
Afin d’illustrer par un exemple les résultats numériques obtenus avec les diverses méthodes, considérons le cas particulier de la fonction et son intégrale sur l’intervalle . Puisque est une primitive de f(x), il est facile de déterminer la valeur exacte de l’intégrale qui vaut I = 1.
En application des formules composites des diverses méthodes, le graphique ci-contre présente le nombre de chiffres significatifs exacts (soit ) obtenus en fonction du nombre n de sous-intervalles de la décomposition.
Remarques :- Le comportement de la fonction choisie étant bien régulier, les diverses courbes croissent très uniformément avec le nombre de sous-intervalles (hormis celles issues de la méthode de Monte-Carlo). Cette apparence est trompeuse : avec une fonction plus chaotique, les courbes seraient beaucoup plus erratiques.
- Pour la méthode de Monte-Carlo, le nombre q indiqué entre parenthèses comme suffixe de la légende (MC (q)) correspond au nombre de tirages réalisés sur chaque sous-intervalle.
- Dans les calculs, la méthode de Romberg n’a pas été traitée dans le cadre où elle présente ses meilleures performances théoriques. Afin de permettre une comparaison avec les autres approches, elle a en effet été appliquée sur chacun des n sous-intervalles utilisés par les autres méthodes[5]. Toutefois, il s’avère que, dans le cas particulier de cet exemple, les deux approches conduisent à des résultats semblables.
Autres méthodes de quadrature numérique
Méthode de calcul d'intégrale à plusieurs dimensions
Méthode de calcul d'intégrale de forme particulière
- Méthode de Laplace pour les intégrales du type ;
- Méthode du point col pour les intégrales du type .
Notes (et références)
- Tronquer l’intervalle pour le rendre fini est une mauvaise idée car la contribution du morceau amputé n’est jamais négligeable.
- Imposer des bornes à la fonction est aussi une mauvaise idée : la contribution élaguée n’est pas négligeable et la régularité est partiellement perdue.
- Les cas singuliers de ce type sont les suivants : (1, 0, r), (1, 2, {1,3,5}), (3, 0, {1,3,5}), (3, 2, {1,3}) et (5, {0,2}, {1,3}).
- compact est uniformément continue et bornée. Toute fonction continue sur un
- n = 2p, on aurait pu appliquer la méthode sur l’intervalle global, mais en effectuant p itérations supplémentaires. Pour le résultat déterminé avec
Voir aussi
Articles connexes
- Intégrale
- Table d'intégrales
- Calcul intégral
- Calcul numérique d'une intégrale
- Primitive
- Règles de Bioche
Catégories :- Algorithme numérique
- Intégration numérique
Wikimedia Foundation. 2010.