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 :

1. Décomposition du domaine en morceaux (un intervalle en sous-intervalles contigus) ;
2. Intégration approchée de la fonction sur chaque morceau ;
3. Sommation des résultats numériques ainsi obtenus.

Définitions

Formule de quadrature

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 Modèle:Math 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.

Ordre d'une méthode de quadrature

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 Modèle:Mvar pour laquelle la valeur approchée de l’intégrale est exacte pour tout polynôme de degré inférieur ou égal à Modèle:Mvar.

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 Modèle:Mvar si la fonction n’est pas continûment dérivable jusqu’à l’ordre Modèle:Math.

Méthode de calcul d'intégrale à une dimension

Généralités

Considérons une intégrale définie <math>I = \int_a^b f(x) \, \mathrm{d}x</math> dont on cherche à estimer la valeur numérique.

Hypothèses et traitements préalables

Supposons que Modèle:Mvar et Modèle:Mvar soient finis : dans le cas contraire, il est conseillé d’effectuer un changement de variable permettant de satisfaire cette hypothèse<ref>Tronquer l’intervalle pour le rendre fini est une mauvaise idée car la contribution du morceau amputé n’est jamais négligeable.</ref>.

Supposons également que la fonction Modèle:Mvar à intégrer comporte une singularité à une des bornes de l'intervalle . Par exemple, la fonction Modèle:Math avec Modèle:Math est intégrable sur Modèle:Math et Modèle:Math. Bien que Modèle:Mvar soit parfaitement régulière sur Modèle:Math, la singularité en 0 et l’impossibilité de la prolonger par une fonction continue causent de grandes difficultés à toutes les méthodes d’intégration numérique, en particulier celles qui utilisent explicitement Modèle:Math<ref>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.</ref>. Dans une telle situation, il convient de soustraire à Modèle:Mvar une fonction Modèle:Mvar dont l’intégrale est connue et qui soit telle que Modèle:Math ne soit plus singulière, puis d’intégrer numériquement cette différence.

Mise en œuvre : décomposition de l'intervalle en morceaux

Considérons une formule de quadrature associée à Modèle:Math du type <math>I_{[0, 1]}(g) = \sum_{i=0}^p \alpha_i \, g(x_i)</math> où les pondérations Modèle:Mvar et les nœuds Modèle:Mvar sont donnés.

Partant d’une décomposition régulière de Modèle:Math en Modèle:Mvar sous-intervalles de longueur Modèle:Math, soit les intervalles Modèle:Math pour Modèle:Math, l’application de la formule de quadrature précédente à chaque Modèle:Mvar s’effectue à l’aide d’une transformation affine, permettant ainsi d’obtenir une approximation Modèle:Math de Modèle:Mvar qui s’écrit :

<math>I_n(f) = h \sum_{k=0}^{n-1} \sum_{i=0}^p \alpha_i \, f(a + h (k + x_i)).</math>

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 Modèle:Math et multiplications par Modèle:Mvar). Cette question est développée plus loin pour quelques formules de quadrature particulières.

Formules utilisant des valeurs des dérivées de la fonction

Si la formule de quadrature comporte des termes du type <math>\sum_{j=0}^{p'} \beta_j \, g'(y_j)</math> faisant intervenir la dérivée ou du type <math>\sum_{j=0}^{p} \gamma_j \, g(z_j)</math> impliquant la dérivée seconde, la transformation affine fait apparaître des facteurs Modèle:Mvar ou Modèle:Math, ceci conformément à la relation suivante :

<math>I_n(f) = h \sum_{k=0}^{n-1} \left( \sum_{i=0}^p \alpha_i \, f(a + h (k + x_i)) + h \sum_{j=0}^{p'} \beta_j \, f'(a + h (k + y_j)) + h^2 \sum_{j=0}^{p} \gamma_j \, f(a + h (k + z_j))\right).</math>

Lien entre ordre de la formule de quadrature et convergence de la méthode

Si Modèle:Mvar est l’ordre de la formule de quadrature et si Modèle:Math est de classe <math>\mathcal{C}_I^{m+1}([a, \, b])</math> (soit l’espace des fonctions Modèle:Math fois dérivables dont la dérivée Modèle:Math est continue par morceaux), notons <math>M = \sup_{x \in [a, b]} |f^{(m+1)}(x)|</math>.

On suppose que la formule de quadrature s'écrit <math>I_n(f) = h \sum_{k=0}^{n-1} \sum_{i=0}^p \alpha_i \, f(a + hk + hx_i)</math>.

Dans ce cas, il existe une constante Modèle:Mvar indépendante de Modèle:Mvar et de Modèle:Math telle que

<math>|I - I_n(f)| \leq C M (b-a) h^{m+1}</math>

Modèle:Démonstration{(m+1)!}</math> pour tout <math>x \in J_k.</math>

Il en découle

<math>|I - I_n(f)| = \left|\sum_{k=0}^{n-1} \int_{J_k} R_{m, k}(x) \, dx - h \sum_{k=0}^{n-1} \sum_{i=0}^p \alpha_i \, R_{m, k}(x_{0, k} + h x_i)\right|</math> <math>\leq \sum_{k=0}^{n-1} \left|\int_{J_k} R_{m, k}(x) \, dx \right| + \left|h \sum_{i=0}^p \alpha_i \, R_{m, k}(x_{0, k} + h x_i)\right|\leq \sum_{k=0}^{n-1} \frac{M h^{m+2}}{(m+1)!} \, \left[\frac{1}{(m+2)} + \sum_{i=0}^p |\alpha_i|\right].</math>

En choisissant <math>\, C = \frac{1}{(m+2)!} \, [1 + (m+2) \sum_{i=0}^p |\alpha_i|],</math> il vient

<math>|I - I_n(f)| \leq C M (b-a) h^{m+1}.</math>

Remarque : ce développement n’a pas pour objectif de déterminer la constante Modèle:Mvar la plus faible. }}

Ce résultat conforte les recommandations suivantes :

• Si Modèle:Math n’est pas continue sur Modèle:Math, une formule de quadrature d’ordre Modèle:Mvar (ou plus) présente peu d’intérêt.
• Si Modèle:Math est régulière par morceaux, il vaut la peine de décomposer Modèle:Math 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

Fichier:Intégration num rectangles.svg

La surface en rouge représente la valeur de l’intégrale estimée par la méthode du point milieu.

C’est la méthode la plus simple qui consiste à interpoler la fonction Modèle:Mvar à intégrer par une fonction constante (polynôme de degré 0).

Si Modèle:Mvar est le point d’interpolation, la formule est la suivante :

<math>I(f) = (b-a) f(\xi)\,</math>

Le choix de Modèle:Mvar influence l’erreur Modèle:Math :

• Si Modèle:Mvar ou Modèle:Mvar, l’erreur est donnée par <math>E(f) = \frac{(b-a)^2}{2} f'(\eta), \quad \eta \in [a,b].</math>C’est la méthode du rectangle qui est d’ordre 0.
• Si Modèle:Math, l’erreur est donnée par<math>E(f) = \frac{(b-a)^3}{24} f(\eta), \quad \eta \in [a,b].</math>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 Modèle:Math) 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 égales et de signes opposés.

Formule du trapèze

Fichier:Intégration num trapèzes.svg

La surface en rouge représente la valeur de l'intégrale estimée par la méthode des trapèzes.

En interpolant Modèle:Mvar par un polynôme de degré 1, les deux points d'interpolation Modèle:Math et Modèle:Math 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 :

<math>I(f) = (b-a) \, \frac{f(a) + f(b)}{2}</math>

L’erreur est donnée par

<math>E(f) = - \frac{(b-a)^3}{12} f(\eta), \quad \eta \in [a,b]</math>

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

Fichier:Intégration num Simpson.svg

La surface en rouge représente la valeur de l'intégrale estimée par la méthode de Simpson.

En interpolant Modèle:Mvar 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 Modèle:Mvar, Modèle:Mvar, et celle choisie en leur milieu Modèle:Math. 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 :

<math>I(f) = \frac{(b-a)}{6} \left[ f(a) + 4 f\left(\frac{a+b}2\right) + f(b) \right]</math>

L’erreur globale est donnée par

<math>E(f) = - \frac{(b-a)^5}{2880} f^{(4)}(\eta), \quad \eta \in [a,b]</math>

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 : formules de Newton-Cotes NC-m

Modèle:Loupe La fonction Modèle:Mvar peut être interpolée à l’aide de son évaluation en Modèle:Mvar points équidistants (comprenant les deux extrémités si Modèle:Math, méthode du point milieu si Modèle:Math) par un polynôme de degré Modèle:Math 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 Modèle:Mvar points, il en découle une méthode

• d’ordre Modèle:Mvar si Modèle:Mvar est impair (« anomalie »),
• d’ordre Modèle:Math si Modèle:Mvar est pair.

On notera ici NC-m la méthode basée sur Modèle:Mvar 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é Modèle:Mvar du polynôme d'interpolation, quitte à subdiviser l'intervalle en sous-intervalles.

Formules 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 Modèle:Mvar points équidistants, mais également la dérivée de la fonction évaluée en Modèle:Mvar points équidistants ; malgré l’abus de langage, on notera ici NC-m-m' une telle formule.

On se limitera ici à Modèle:Math correspondant aux deux extrémités Modèle:Mvar et Modèle:Mvar.

Formules NC-m-2

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 Modèle:Math et de Modèle:Math 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 Modèle:Mvar désigne les points d’évaluation de Modèle:Mvar (Modèle:Mvar entre 0 et Modèle:Math) :

• Formule NC-1-2 : basée sur un polynôme de degré 2, elle est d’ordre 3 :

<math>I(f) = \frac{(b-a)}{24} \left[ 24 f \left(\frac{b+a}{2}\right) +(b-a)(f'(b)-f'(a)) \right]</math>

• Formule NC-2-2 : basée sur un polynôme de degré 3, elle est d’ordre 3 :

<math>I(f) = \frac{(b-a)}{12} \left[ 6 (f(a) + f(b)) + (b-a)(f'(a) - f'(b)) \right]</math>

• Formule NC-3-2 : basée sur un polynôme de degré 4, elle est d’ordre 5 :

<math>I(f) = \frac{(b-a)}{60} \left[ 14(f(a) + f(b)) + 32 f(x_1) + (b-a)(f'(a) - f'(b)) \right]</math>

• Formule NC-4-2 : basée sur un polynôme de degré 5, elle est d’ordre 5 :

<math>I(f) = \frac{(b-a)}{240} \left[ 39 (f(a) + f(b)) + 81 (f(x_1) + f(x_2)) + (b-a)(f'(a) - f'(b)) \right]</math>

• Formule NC-5-2 : basée sur un polynôme de degré 6, elle est d’ordre 7 :

<math>I(f) = \frac{(b-a)}{3780} \left[ 434 (f(a) + f(b)) + 1024 (f(x_1) + f(x_3)) + 864 f(x_2) + (b-a)(f'(a) - f'(b)) \right]</math>

Concernant l’erreur globale d’une formule de quadrature linéaire d’ordre Modèle:Mvar, elle est donnée par

<math>E(f) = C (b-a)^{p+2} f^{(p+1)}(\eta), \quad \eta \in [a,b]</math>

Ce procédé permet ainsi une généralisation des formules de Newton-Cotes. Avec Modèle:Mvar points pour la fonction et 2 points pour sa dérivée, il en découle une méthode

• d’ordre Modèle:Math si Modèle:Mvar est impair (anomalie),
• d’ordre Modèle:Math si Modèle:Mvar est pair.

Formules NC-m-m'-m"

La notation indique que la dérivée seconde intervient également en Modèle:Math points équidistants. Mentionnons uniquement une formule particulièrement remarquable qui présente une double anomalie :

Formule NC-4-2-2, basée sur un polynôme de degré 7, elle est d’ordre 9 :

<math>I(f) = \frac{(b-a)}{6720} [ 1173 (f(a) + f(b)) + 2187 (f(x_1) + f(x_2)) + 78 (b-a)(f'(a) - f'(b)) + 2 (b-a)^2(f(a) + f(b)) ]</math>

Les coefficients de Modèle:Math et de Modèle:Math sont opposés, ce qui permet des simplifications dans la formule composite. Par contre, ce n’est pas le cas pour Modèle:Mvar.

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'-mModèle:'', m doit être non nul afin de tenir compte d’une translation sur Modèle:Math.
• Tout triplet d’entiers (m, m', mModèle:'') ne conduit pas nécessairement à une formule de quadrature. C’est le cas pour les triplets du type (2p, 2q +1, r), ainsi que quelques cas du type (2p + 1, 2q, 2r +1) <ref>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}).</ref>.
• Les triplets du type (2p +1, 2q, 2 r), (2p +1, 2q + 1, 2r + 1) et (2p, 2q, 2r +1) présentent une anomalie dans le sens où la formule est basée sur un polynôme de degré Modèle:Math (nombre de degrés de liberté) alors qu’elle est encore vraie pour un polynôme de degré Modèle:Math.

Formules composites

Pour chacune des méthodes précédentes, le terme d’erreur dépend d’une puissance de Modèle:Math. Cette imprécision étant le plus souvent trop importante, l’erreur peut être réduite en découpant simplement l’intervalle Modèle:Math 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 Modèle:Math est estimée par la somme des valeurs ainsi calculées.

On appelle formule composite l’expression caractérisant cette estimation.

Notons :

• Modèle:Mvar l’indice des Modèle:Mvar sous-intervalles ;
• Modèle:Math la longueur de chacun d’eux ;
• Modèle:Math la borne inférieure ;
• Modèle:Math le point milieu,

ceci pour Modèle:Mvar entre 0 et Modèle:Math. Les formules composites sont les suivantes :

• Méthode du point milieu d’ordre 1 :

<math> I(f) = \frac{b-a}{n} \sum_{k=0}^{n-1} f( m_k )</math>

• Méthode des trapèzes d’ordre 1 :

<math>I(f) = \frac{b-a}{n} \left( {f(a) + f(b) \over 2} + \sum_{k=1}^{n-1} f ( a + k h) \right)</math>

• Méthode de Simpson d’ordre 3 :

<math>I(f) = \frac{h}{6} \left( f(a) + f(b) + 2 \sum_{k=1}^{n-1} f(x_k) + 4 \sum_{k=0}^{n-1} f(m_k) \right)</math>

Pour les formules de quadrature faisant intervenir des dérivées, voici deux exemples, les autres cas se déduisant par analogie :

• NC-3-2 d’ordre 5 :

<math>I(f) = \frac{h}{60} \left( 14 (f(a) + f(b)) + h (f'(a) - f'(b)) + 28 \sum_{k=1}^{n-1} f(x_{k}) + 32 \sum_{k=0}^{n-1} f(m_k) \right)</math>

• NC-4-2-2 d’ordre 9 :

<math>I(f) = \frac{h}{6720} \left( 1173 (f(a) + f(b))+ 78 h (f'(a) - f'(b)) + 2 h^2 (f(a) + f(b))\right.</math> <math> \left. + 2346 \sum_{k=1}^{n-1} f(x_{k}) + 2187 \sum_{k=0}^{n-1} (f(x_{k, 1}) + f(x_{k, 2})) + 4 h^2 \sum_{k=1}^{n-1} f(x_{k}) \right)</math>

(Ici Modèle:Math où Modèle:Math, Modèle:Math étant le nombre de points d’un sous-intervalle)

Concernant l’erreur finale d’une formule de quadrature linéaire d’ordre Modèle:Mvar, elle est donnée par la relation

<math>E_h(f) = C h^{p+1} (b-a) f^{(p+1)}(\eta),\quad \eta\in [a,b]</math>

Méthode de Romberg

Modèle:Article détaillé 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 Modèle:Mvar itérations, elle conduit à une méthode d’ordre de Modèle:Math avec une erreur en Modèle:Math pour des fonctions de classe <math>\mathcal{C}_I^{2n+2}</math>, ceci avec une grande économie du nombre d’évaluations de la fonction (précisément Modèle:Math évaluations).

Tirages aléatoires

Méthode de Monte-Carlo

Modèle:Article détaillé Pour intégrer une fonction Modèle:Mvar sur un intervalle Modèle:Math, 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 Modèle:Mvar qui sont nécessaires pour espérer obtenir un résultat significatif.

L’estimation de l’intégrale Modèle:Mvar de Modèle:Mvar est fournie par Modèle:Math où Modèle:Math est la moyenne arithmétique des Modèle:Math, les Modèle:Mvar étant Modèle:Mvar tirages aléatoires indépendants et distribués uniformément sur Modèle:Math.

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 Modèle:Mvar tirages correspond une estimation particulière de l’intégrale, c'est-à-dire une réalisation d’une variable aléatoire Modèle:Math dont la distribution dépend de Modèle:Mvar, de Modèle:Math et de Modèle:Mvar. On peut alors assurer que son espérance est égale à l’intégrale de Modèle:Mvar et préciser une borne de l’écart type de l’erreur.

Distribution de l'erreur

Les résultats suivants permettent de caractériser la distribution de l’erreur et son écart type :

• Lorsque Modèle:Mvar augmente, l’erreur Modèle:Math s’amenuise statistiquement. Plus précisément :

Soit Modèle:Mvar l’écart type d’un tirage qui est défini par la relation

<math>\sigma^2 = \frac{1}{b - a} \int_a^b f^2(x) \,\mathrm{d}x - \left( \frac{I}{b - a} \right)^2 .</math>

Alors la distribution de <math display=inline>\frac {\sqrt{q}}{(b - a) \sigma} \, E_q(f)</math> converge vers celle de la loi normale centrée réduite, soit <math>\mathcal{N}(0,\, 1)</math>. En d’autres termes, Modèle:Math est essentiellement distribué comme une loi normale <math display=inline>\mathcal{N}\left(0,\, \frac {(b - a)^2}{q} \sigma^2 \right).</math>

• Si Modèle:Mvar est une fonction bornée quelconque (la borne est notée <math>\|f\|_\infty</math>), alors <math>\sigma \leq \|f\|_\infty</math> et ainsi

<math>\sigma (E_q(f)) \leq \frac {(b-a)}{\sqrt{q}} \, \|f\|_\infty</math>.

• Si Modèle:Mvar est de classe <math>\mathcal{C}^1([a, \, b])</math> (dérivable et à dérivée continue, en particulier bornée<ref>Toute fonction continue sur un compact est uniformément continue et bornée.</ref>), alors <math display=inline>\sigma \leq \frac{b - a}{\sqrt{3}} \|f'\|_\infty</math> et ainsi

<math>\sigma (E_q(f)) \leq \frac {(b-a)^2}{\sqrt{3 q}} \, \|f'\|_\infty.</math>

• Si Modèle:Mvar est de classe <math>\mathcal{C}^1([a, \, b])</math>, que Modèle:Math est décomposé en Modèle:Mvar intervalles de longueur Modèle:Mvar avec Modèle:Math, que le nombre total Modèle:Mvar de tirages est réparti à parts égales dans chaque intervalle, alors

<math>\sigma (E_q(f)) \leq (b-a) \, \frac {h}\sqrtModèle:3 q \, \|f'\|_\infty.</math>

Modèle:Démonstration \|f'\|_\infty.</math>

En utilisant l’indépendance des variables Modèle:Math et l’hypothèse que les Modèle:Mvar sont tous égaux à Modèle:Math, il vient

<math>\sigma^2 (E_q(f)) = \sum_{i = 1}^n \sigma^2 (E_{J_i, q_i}(f)) \leq (b - a)^2 \frac{h^2}{3 q} \|f'\|_\infty^2.</math>

}}

Ce dernier résultat justifie le comportement relativement médiocre de la méthode de Monte-Carlo. Supposons en effet que le nombre Modèle:Mvar de tirages soit imposé (limitation de l’effort de traitement). Dans ce cas et afin de réduire <math>\sigma (E_q(f))</math>, il convient de choisir le plus petit Modèle:Mvar possible, soit un tirage par intervalle Modèle:Math pour obtenir

<math>\sigma (E_q(f)) \leq \sqrt {\frac{(b-a)}{3}} h^{3/2} \, \|f'\|_\infty.</math>

La méthode du point milieu est plus efficace si la fonction est dans <math>\mathcal{C}^2([a, \, b])</math>.

Tableau comparatif

Le tableau suivant résume les performances théoriques de chaque méthode :

Nom de la méthode	Degré du polynôme	Nombre de points	Degré d’exactitude (ordre)	Degré d’erreur globale	Degré d’erreur finale
Rectangle	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
Simpson 3/8	3	4	4	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 Modèle:Mvar dans l’erreur pour une application globale de la formule sur l’intervalle.
• Degré d’erreur finale : la puissance du facteur Modèle:Mvar dans la formule composite.

Exemple numérique

Fichier:Integration num.png

Performances comparatives de quelques méthodes d’intégration : résultats tirés d’un exemple.

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 Modèle:Math et son intégrale sur l’intervalle <math>[0, (\pi)^{\frac{1}{2}}]</math>. Puisque <math>F(x) = - \frac{1}{2} \, \cos(x^2)</math> est une primitive de Modèle:Math, il est facile de déterminer la valeur exacte de l’intégrale qui vaut Modèle:Math.

En application des formules composites des diverses méthodes, le graphique ci-contre présente le nombre de chiffres significatifs exacts (soit <math>- \log_{10}\left(\left|\frac{E(f)}{I}\right|\right)</math>) obtenus en fonction du nombre Modèle:Mvar 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 Modèle:Mvar 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 Modèle:Mvar sous-intervalles utilisés par les autres méthodes<ref>Pour le résultat déterminé avec Modèle:Math, on aurait pu appliquer la méthode sur l’intervalle global, mais en effectuant Modèle:Mvar itérations supplémentaires.</ref>. Toutefois, il s’avère que, dans le cas particulier de cet exemple, les deux approches conduisent à des résultats semblables.

Méthodes de Gauss

Dans le cas du calcul d'une intégrale de la forme <math>\int_a^b f(x) \omega(x) \,\mathrm{d}x</math>, où Modèle:Math est vue comme une fonction poids déterminée sur Modèle:Math (les bornes peuvent être infinies), les méthodes de quadrature de Gauss permettent de calculer efficacement cette intégrale en utilisant les propriétés des polynômes orthogonaux pour Modèle:Math.

Méthode de calcul d'intégrale de forme particulière

• Méthode de Laplace pour les intégrales du type <math>\int_a^b\! e^{M f(x)}\, \mbox{d}x\,</math>;
• Méthode du point col pour les intégrales du type <math>I(\lambda) = \int_\mathcal{C} f(z) e^{\lambda g(z)} \, \mbox{d}z\,</math>.

Calculs d'intégrales sur intervalles infinis

On peut adapter le principe des méthodes de calcul numérique d'intégrale aux cas où l'intervalle d'intégration n'est pas borné. On peut appliquer des méthodes de quadrature adaptées, comme la méthode de quadrature de Gauss-Hermite pour une intégrale sur la droite réelle ou de Gauss-Laguerre sur la demi-droite réelle positive<ref>Modèle:Ouvrage</ref>. On peut aussi envisager des méthodes de Monte-Carlo, ou une transformation qui permet de revenir à un intervalle fermé, par exemple : <math display="block">

\int_{-\infty}^{\infty} f(x) \,\mathrm{d}x

= \int_{-1}^{+1} f\left( \frac{t}{1-t^2} \right) \frac{1+t^2}{\left(1-t^2\right)^2} \,\mathrm{d}t, </math> ou pour le cas semi-infini <math display="block">\begin{align} \int_a^{\infty} f(x) \,\mathrm{d}x &= \int_0^1 f\left(a + \frac{t}{1-t}\right) \frac{\mathrm{d}t}{(1-t)^2}, \\ \int_{-\infty}^a f(x) \,\mathrm{d}x &= \int_0^1 f\left(a - \frac{1-t}{t}\right) \frac{\mathrm{d}t}{t^2}, \end{align}.</math>

Méthode de calcul d'intégrale à plusieurs dimensions

• Méthode de Monte-Carlo

Évaluation de l'erreur

Noyau de Peano d'une méthode

En se plaçant dans le cadre général d'une méthode d'intégration numérique, l'erreur commise s'exprime sous la forme :

<math>E(f) = \int_a^b f(t) w(t) \,\mathrm{d}t - \sum_{i=1}^n \omega_i f(\xi_i).</math>

L'erreur peut être calculée par l'intégrale d'une fonction liée à l'erreur, appelée noyau de Peano : Modèle:Théorème

L'intégrale est une forme linéaire en la fonction Modèle:Mvar. Le noyau de Peano permet de contrôler l'erreur, notamment par les deux résultats suivants : Modèle:Énoncé Modèle:Énoncé

Exemples

On donne ici quelques calculs du noyau de Peano.

Noyau de la méthode du point médian

L'erreur pour cette méthode d'ordre 1 s'écrit :

<math>K_1(t) = \int_a^b (x-t)_+ \,\mathrm{d}x - (b-a) \left(\frac{a+b}{2}-t\right)_+ = \int_t^b (x-t) \,\mathrm{d}x - (b-a) \left(\frac{a+b}{2}-t\right)_+</math>

<math>K_1(t)= \frac12 (b-t)^2 - (b-a) \left(\frac{a+b}{2}-t\right)_+.

</math>

<math>K_1(t) = \begin{cases} \frac12 (t-a)^2 & \text{ si } t \leq \tfrac{a+b}{2}, \\ \frac12 (b-t)^2 & \text{ si } t \geq \tfrac{a+b}{2}.\end{cases}

</math> Il vient alors, par application du deuxième résultat :

<math>E(f) = \frac{1}{1!} f^{}(\xi) \int_a^b K_1(t) \,\mathrm{d}t = \frac{(b-a)^3}{24}f^{}(\xi) .</math>

Noyau de la méthode des trapèzes

Le noyau de Peano pour cette méthode d'ordre 1 s'écrit :

<math>K_1(t) = \frac12 (t-a)(t-b),</math>

Erreur de la méthode de quadrature de Gauss

On sait que la méthode de quadrature de Gauss de degré m est d'ordre 2m+1. Il n'existe pas de formule générale dans ce cas, mais on peut obtenir le résultat suivant<ref>Modèle:Ouvrage</ref> : Modèle:Énoncé

Articles connexes

QSS, méthode d'intégration à événements discrets