Mat 367, Méthodes numériques

Bernard.Parisse@ujf-grenoble.fr

2016

La version HTML de ce document comporte des champs de saisie interactifs, ceux-ci apparaissent comme des commandes “mortes” dans la version PDF. La version HTML est optimisée pour le navigateur Firefox. Vous pouvez exécuter toutes les commandes interactives en cliquant sur le bouton Exécuter (attention cela peut prendre un certain temps!), le champ suivant est la console de l’interpréteur du logiciel de calcul formel.
Les commandes sont exécutées par la version javascript de Giac/Xcas. Bien que le thème de ce cours soit numérique, on verra sur de nombreuses commandes qu’il est fort utile de disposer d’un logiciel capable de faire du calcul formel!

Table des matières

Index

1  Présentation du module

Les thèmes abordés seront :

  1. calcul approché, représentation des données (flottants, vecteurs, matrices), erreurs (normes, de calcul, d’arrondi...).
  2. Pivot de Gauss, factorisation LU, conditionnement, Cholesky, factorisation QR.
  3. Interpolation polynômiale (évaluation, interpolation de Lagrange, Hermite, Bézier)
  4. Intégration numérique
  5. Méthode du point fixe, de Newton, méthodes itératives en algèbre linéaire
  6. Méthode de la puissance, valeurs propres et vecteurs propres.
  7. Résolution d’équations différentielles.

L’évaluation se fait sur :

Les calculatrices et les netbooks de taille d’écran plus petits que 13 pouces sont autorisées au DS et à l’examen final (prêt possible de netbooks pour le semestre).

2  Représentation des nombres et autres données, calcul exact/approché

Mot-clefs:
Types de base : entier machine, entier long, flottant machine et multiprécision (Base 2, base 10).
Erreur relative, erreur absolue, erreur d’arrondi, +/-, */%
Algorithme de Horner. Types composés : complexes, polynomes (représentation dense/creuse), symboles, listes (vecteurs, matrices), expressions, fonctions.
Les principaux ensembles de nombres en mathématiques sont les entiers positifs ℕ et relatifs ℤ, les rationnels ℚ, les réels ℝ et les complexes ℂ. Sur ordinateur, on peut représenter ces nombres de manière exacte dans certains cas, approchée dans d’autres.

2.1  Représentation des entiers

Proposition 1   Division euclidienne de deux entiers : si a et b sont deux entiers, a ≥ 0, b>0, il existe un unique couple (q,r) tel que
a = bq +r ,    r ∈ [0, b

Preuve : On prend pour q le plus grand entier tel que abq ≥ 0.
Exemple :

La division euclidienne permet d’écrire un nombre entier, en utilisant une base b et des caractères pour représenter les entiers entre 0 et b−1. Nous écrivons les nombres entiers en base b=10 avec comme caractères les chiffres de 0 à 9. Les ordinateurs utilisent des circuits binaires pour stocker les informations, il est donc naturel d’y travailler en base 2 en utilisant comme caractères 0 et 1 ou en base 16 en utilisant comme caractères les chiffres de 0 à 9 et les lettres de A à F. En général, pour trouver l’écriture d’un nombre en base b (par exemple b=2), on effectue des divisions euclidienne successives par b du nombre puis de ses quotients successifs jusqu’à ce que le quotient soit 0 et on accolle les restes obtenus (premier reste à droite, dernier reste à gauche). Inversement, pour retrouver un entier d à partir de son écriture dn...d0, on traduit les divisions euclidiennes successives en

 d=( ... ((dn b +dn−1)b + dn−2)...+d1)b+d0
 =dn bn + dn−1 bn−1 + ... + d0

Par exemple, vingt-cinq s’écrit en base 16 0x19 car 25 divisé par 16 donne quotient 1, reste 9


En base 2, on trouverait 0b11001 car 25=24+23+1.


On peut effectuer les opérations arithmétiques de base (+,-,*, division) directement en base 2 (ou 16). Par exemple la table de l’addition est 0+0=0, 0+1=1+0=1 et 1+1=0 je retiens 1, donc :

  01001111
+ 01101011
----------
  10111010

Exercice : comment passe-t-on simplement de la représentation d’un nombre en base 2 à un nombre en base 16 et réciproquement ?

Les microprocesseurs peuvent effectuer directement les opérations arithmétiques de base sur les entiers “machine” (déclinés en plusieurs variantes selon la taille et la possibilité d’avoir un signe). Noter que la division de deux entiers a et b n’a pas la même signification que la division de deux réels, comme elle ne tomberait pas forcément juste, on calcule le quotient et le reste de la division euclidienne.

Ces entiers machines permettent de représenter de manière exacte des petits entiers relatifs par exemple un entier machine signé sur 4 octets est compris entre [−231,231−1].

Ces entiers machines permettent de faire très rapidement du calcul exact sur les entiers, mais à condition qu’il n’y ait pas de dépassement de capacité, par exemple pour des entiers 32 bits, 230+230+230+230 renverra 0. Ils sont utilisables avec tous les langages de programmation traditionnels.

Les logiciels de calcul formel et certains logiciels de programmation permettent de travailler avec des entiers de taille beaucoup plus grande, ainsi qu’avec des rationnels, permettant de faire du calcul exact, mais on paie cette exactitude par un temps de calcul plus long, de plus pas mal de méthodes numériques ne gagnent rien à faire des calculs intermédiaires exacts. Néanmoins, l’utilisation d’un logiciel de calcul formel permettra dans certains cas d’illustrer certains phénomènes dus au calcul approché.

2.2  Les réels

On se ramène d’abord au cas des réels positifs, en machine on garde traditionnellement un bit pour stocker le signe du réel à représenter.

2.2.1  Virgule fixe et flottante.

La première idée qui vient naturellement serait d’utiliser un entier et de déplacer la virgule d’un nombre fixe de position, ce qui revient à mulitplier par une puissance (négative) de la base. Par exemple en base 10 avec un décalage de 4, 1234.5678 serait représenté par 12345678 et 1.2345678 par 12345 (on passe de l’entier au réel par multiplication par 10−4. L’inconvénient d’une telle représentation est qu’on ne peut pas représenter des réels grands ou petits, comme par exemple le nombre d’Avogadro, la constante de Planck, etc.

D’où l’idée de ne pas fixer la position de la virgule, on parle alors de représentation à virgule flottante ou de nombre flottant : on représente un nombre par deux entier, l’un appelé mantisse reprend les chiffres significatifs du réel sans virgule, l’autre l’exposant, donne la position de la virgule. Attention, le séparateur est un point et non une virgule dans la grande majorité des logiciels scientifiques. On sépare traditionnellement la mantisse de l’exposant par la lettre e. Par exemple 1234.5678 peut être représenté par 12345678e-8 (mantisse 12345678, exposant -8) mais aussi par 1234567800e-10.

Naturellement, sur un ordinateur, il y a des limites pour les entiers représentant la mantisse m et l’exposant e. Si on écrit les nombres en base b, la mantisse m s’écrira avec un nombre n fixé de chiffres (ou de bits en base 2), donc m ∈ [0,bn[. Soit un réel x représenté par

x=mbe,    m ∈ [0,bn

Si m∈ [0,bn−1[, alors on peut aussi écrire x=mbe−1 avec m′=mb ∈ [0,bn[, quelle écriture faut-il choisir? Intuitivement, on sent qu’il vaut mieux prendre m′ le plus grand possible, car cela augmente le nombre de chiffres significatifs (alors que des 0 au début de m ne sont pas significatifs). Ceci est confirmé par le calcul de l’erreur d’arrondi pour représenter un réel. En effet, si x est un réel non nul, il ne s’écrit pas forcément sous la forme mbe, on doit l’arrondir, par exemple au plus proche réel de la forme mbe. La distance de x à ce réel est inférieure ou égale à la moitié de la distance entre deux flottants consécutifs, mbe et (m+1)be, donc l’erreur d’arrondi est inférieure ou égale à be/2. Si on divise par xmbe, on obtient une erreur relative d’arrondi majorée par 1/(2m). On a donc intérêt à prendre m le plus grand possible pour minimiser cette erreur. Quitte à mulitplier par b, on peut toujours se ramener (sauf exceptions, cf. ci-dessous), à m ∈ [bn−1,bn[, on a alors une erreur d’arrondi relative majorée par

1
2bn−1

On appelle flottant normalisé un flottant tel que m ∈ [bn−1,bn[. Pour écrire un réel sous forme de flottant normalisé, on écrit le réel en base b, et on déplace la virgule pour avoir exactement n chiffres non nuls avant la virgule et on arrondit (par exemple au plus proche). L’exposant est égal au décalage effectué. Notez qu’en base 2, un flottant normalisé commence forcément par 1, ce qui permet d’économiser un bit dans le stockage.

Ainsi, l’erreur d’arrondi commise lorsqu’on représente un réel (connu exactement) par un double normalisé est une erreur relative inférieure à de 2−53 (b=2 et n=52+1 pour les doubles).

Exemples :

Il existe une exception à la possibilité de normaliser les flottants, lorsqu’on atteint la limite inférieure de l’exposant e. Soit en effet em le plus petit exposant des flottants normalisés et considérons les flottants x=bem(1+1/b) et y=bem. Ces flottants sont distincts, mais leur différence n’est plus représentable par un flottant normalisé. Comme on ne souhaite pas représenter xy par 0, (puisque le test x==y renvoie faux), on introduit les flottants dénormalisés , il s’agit de flottants dont l’exposant est l’exposant minimal représentable sur machine et dont la mantisse appartient à [0,bn−1[. Par exemple 0 est représenté par un flottant dénormalisé de mantisse 0 (en fait 0 a deux reprśentation, une de signe positif et une de signe négatif).

Enfin, on utilise traditionnellement une valeur de l’exposant pour représenter les nombres plus grands que le plus grand réel reprśentable sur machine (traditionnellement appelé plus ou moins infini) et les erreurs (par exemple 0./0. ou racine carrée d’un nombre réel négatif, traditionnellement appelé NaN, Not a Number).

Exercice : quels sont les nombres réels représentables exactement en base 10 mais pas en base 2 ? Si on écrit 1/10 en base 2 avec 53 bits de précision, puis que l’on arrondit avec 64 bits de précision, ou si on écrit 1/10 en base 2 avec 64 bits de précision, obtient-on la même chose ?

Les ordinateurs reprśentent généralement les flottants en base 2 (cf. la section suivante pour plus de précisions), mais cette représentation n’est pas utilisée habituellement par les humains, qui préfèrent compter en base 10. Les ordinateurs effectuent donc la conversion dans les routines d’entrée-sortie. Le format standard utilisé pour saisir ou afficher un nombre flottant dans un logiciel scientifique est composé d’un nombre à virgule flottante utilisant le point comme séparateur décimal (et non la virgule) suivi si nécessaire de la lettre e puis de l’exposant, par exemple 1.23e-5 ou 0.0000123. Dans les logiciels de calcul formel, pour distinguer un entiers représentés par un entier d’un entier représenté par un flottant on écrit l’entier suivi de .0 par exemple 23.0.

Remarque :
Les microprocesseurs ayant un mode BCD peuvent avoir un format de représentation des flottants en base 10, les nombres décimaux comme par exemple 0.3 peuvent être représentés exactement. Certains logiciels, notamment maple, utilisent par défaut des flottants logiciels en base 10 sur des microprocesseurs sans mode BCD, ce qui entraine une baisse de rapidité importante pour les calculs numériques (on peut partiellement améliorer les performances en utilisant evalhf en maple).

2.2.2  Les flottants au format double

Cette section développe les notions de la section précédente pour les flottants machine selon la norme IEEE-754, utilisables dans les langage de programmation usuels, elle peut être omise en première lecture. La représentation d’un double en mémoire se compose de 3 parties : le bit de signe s=± 1 sur 1 bit, la mantisse M ∈ [0,252[ sur 52 bits, et l’exposant e ∈ [0, 211[ sur 11 bits. Pour les nombres “normaux”, l’exposant est en fait compris entre 1 et 211−2, le nombre représenté est le rationnel

(1+
M
252
) 2e+1−210 

Pour écrire un nombre sous cette forme, il faut d’abord chercher par quel multiple de 2 il faut le diviser pour obtenir un réel r dans [1,2[, ce qui permet de déterminer l’exposant e. Ensuite on écrit la représentation en base 2 de r−1 ∈ [0,1[. Exemples :

On observe que la représentation en base 2 de 6.4 a du être arrondie (car elle est infinie en base 2) bien qu’elle soit exacte (finie) en base 10. Seuls les entiers et les rationnels dont le dénominateur est une puissance de 2 peuvent être représentés exactement. Ceci entraine des résultats qui peuvent surprendre comme par exemple le fait que 0.5 - 5*0.1 n’est pas nul.

Des représentations spéciales (avec e=0 ou e=211−1) ont été introduites pour représenter ± ∞ (pour les flottants plus grands en valeur absolue que le plus grand flottant représentable), et pour représenter les nombres non nuls plus petits que le plus petit flottant représentable de la manière exposée ci-dessus (on parle de flottants dénormalisés), ainsi que le nombre NaN (Not a Number) lorsqu’une opération a un résultat indéfini (par exemple 0/0).

Remarque : Sur les processeurs compatibles avec les i386, le coprocesseur arithmétique i387 gère en interne des flottants avec 80 bits dont 64 bits de mantisse. Sur les architectures 64 bits (x86 ou AMD), le jeu d’instruction SSE permet de travailler avec des flottants de 128 bits. Le compilateur gcc permet d’utiliser ces flottants longs avec le type long double ou les types __float80 et __float128 en utilisant un drapeau de compilation du type -msse

2.2.3  Opérations sur les flottants

Les opérations arithmétiques de base sur les flottants se font de la manière suivante :

2.2.4  Erreurs

La représentation des nombres réels par des doubles présente des avantages, les opérations arithmétiques sont faites au plus vite par le microprocesseur. Les coprocesseurs arithmétiques (intégrés sur les microprocesseurs de PC) proposent même le calcul des fonctions usuelles (trigonométriques, racine carrée, log et exp) sur le type double et utilisent des formats de représentation interne ayant plus de 64 bits pour les doubles, ce qui permet de limiter les erreurs d’arrondi. Par contre, des erreurs vont être introduites, on parle de calcul approché par opposition au calcul exact sur les rationnels. En effet, la représentation doit d’abord arrondir tout réel qui n’est pas un rationnel dont le dénominateur est une puissance de 2. Ensuite chaque opération va entrainer une propagation de ces erreurs et va y ajouter une erreur d’arrondi sur le résultat. Enfin, l’utilisation du type double peut provoquer un dépassement de capacité (par exemple 100!*100!).

Pour diminuer ces erreurs et les risques de dépassement de capacité, il existe des types flottants multiple précision, qui permettent de travailler avec un nombre fixé à l’avance de décimales et une plage d’exposants plus grande. Les calculs sont plus longs mais les erreurs plus faibles. Attention, il s’agit toujours de calcul approché! De plus, pour des quantités dont la valeur est déterminée de manière expérimentale, la source principale de propagation d’erreurs est la précision des quantités initiales, il ne sert souvent à rien d’utiliser des types flottants multiprécision car les erreurs dus à la représentation (double) sont négligeables devant les erreurs de mesure. Dans ce cas, il est pertinent lorsqu’on évalue f(x) avec x mal connu de calculer aussi f′(x), en effet :

f(x(1+h))= f(x)+xh f′(x) + O(h2)

l’erreur relative sur f(x) est donc au premier ordre multipliée par

|
xf′(x)
f(x)
|

Par exemple, l’erreur relative sur ex est au premier ordre l’erreur relative sur x multipliée par |x|.


2.2.5  Erreur absolue, relative, arrondi propagation des erreurs.

On a vu précédemment que pour représenter un réel, on devait l’arrondir, ce qui introduit une erreur même si le réel est connu exactement (par exemple 1/10). Voyons comment se propagent les erreurs dans les opérations arithmétiques de base : on distingue l’addition, la multiplication et l’inversion. La soustraction se ramène à l’addition car le calcul de l’opposé n’introduit aucune erreur nouvelle. Pour l’addition, si |xx0| ≤ ε0 et si |yy0| ≤ ε1 alors par l’inégalité triangulaire (|a+b|≤ |a|+|b|), on a :

|(x+y)−(x0+y0)| ≤ |xx0| + | yy0 | ≤  ε0 + ε1 

on dit que les erreurs absolues s’additionnent.

Définition 1   L’erreur absolue est définie comme un majorant de la valeur absolue de la différence entre le nombre réel et son représentant double :
|xx0| ≤ ε 

Mais comme il faut représenter x0+y0 en machine, on doit ajouter une erreur d’arrondi, qui est proportionnelle à la valeur absolue de x0+y0 d’où la notion d’erreur relative :

Définition 2   L’erreur relative est égale à l’erreur absolue divisée par la valeur absolue du nombre
|xx0| ≤ ε |x0

Remarquons au passage que les erreurs de mesure expérimentales sont pratiquement toujours des erreurs relatives.

Donc lorsqu’on effectue une addition (ou une soustraction) de deux réels sur machine, on doit additionner les deux erreurs absolues sur les opérandes et ajouter une erreur d’arrondi (relative de 2−53, à titre d’exercice, on pourra vérifier que cette erreur d’arrondi est majorée par l’erreur absolue de la somme x+y dès l’instant où x et y ont eux-même une erreur d’arrondi).

Lorsqu’on effectue une multiplication de deux nombres x,y dont les représentants x0,y0 sont non nuls, on a




xyx0 y0
x0 y0
 





x
x0
 
y
y0
 −1 





(
x
x0
−1)(
y
y0
 −1)+(
x
x0
−1)+(
y
y0
 −1) 


l’erreur relative est donc la somme des erreurs relatives et du produit des erreurs relatives (on peut souvent négliger le produit devant la somme). Il faut aussi y ajouter une erreur relative d’arrondi de 2−53 sur x0 y0.

On observe que la multiplication est une opération posant moins de problèmes que l’addition, car on manipule toujours des erreurs relatives, par exemple si l’erreur relative sur deux doubles x et y non nuls est de 2−53, alors l’erreur relative sur xy sera de

2−53 + 2−53 + 2−106 + 2−53 ≈ 3 × 2−53 

Lorsque l’erreur relative sur les données est grande devant 2−53, l’erreur relative d’arrondi final est négligeable, on peut alors dire que les erreurs relatives s’additionnent pour un produit (c’est aussi vrai pour un quotient: exercice!). Par contre, si on additionne deux nombres dont le représentant de la somme est proche de 0, la somme des erreurs absolues peut devenir non négligeable par rapport à la somme des représentants, entrainant une erreur relative très grande. Par exemple si x est représenté par x0=1+2−52 avec une erreur d’arrondi de 2−53 et y par y0=−1 avec la même erreur d’arrondi, l’addition de x et y renvoie 2−52 avec une erreur absolue de 2 * 2−53 (ici il n’y a pas d’arrondi lorsqu’on fait la somme). C’est une erreur relative de 1 (qui domine largement l’erreur d’arrondi) ce qui signifie que dans la mantisse, seul le premier bit sur les 52 a un sens, la perte de précision est très grande.

Une autre conséquence importante est que l’addition de réels sur machine n’est pas une opération associative, par exemple

(2.0−53+2.0−53)+1.0 → 1+2−52 

alors que

2.0−53+(2.0−53+1.0) → 1 

Dans Xcas, il n’y a que 48 bits de mantisse :

Si on a plusieurs termes à additionner, il faut commencer par additionner entre eux les termes les plus petits, pour que les petits termes ne soient pas absorbés un à un dans les erreurs d’arrondi (les petits ruisseaux font les grands fleuves).

Exercice : pour calculer la valeur numérique d’une dérivée de fonction, il vaut mieux calculer (f(x+h)−f(xh))/(2h) que (f(x+h)−f(x))/h car le terme d’erreur est en O(h2) et non en O(h). Attention toutefois à ne pas prendre h trop petit, sinon x+h=x en flottants et même si x+hx, l’erreur absolue sur f(x+h)−f(xh) est (au moins) d’ordre ε |f(x)|, donc l’erreur relative est d’ordre ε/h |f(x)|. Par exemple pour h=1e-8 le reste est en O(h2) donc de l’ordre des erreurs d’arrondi mais l’erreur relative sur f(x+h)−f(xh) est d’ordre є/h largement supérieure (en flottants double-précision). On choisira plutôt h tel que є/h soit proche de h2, donc de l’ordre de 1e-5, qui fournira une valeur approchée avec une erreur relative de l’ordre de 1e-10. Exemple : calcul de la dérivée numérique de exp(sin(x)) en x=1

Remarquons néanmoins que les erreurs calculées ici sont des majorations des erreurs réelles (ou si on préfère l’erreur obtenue dans le pire des cas), statistiquement les erreurs sur les résultats sont moindres, par exemple si on effectue n calculs susceptibles de provoquer des erreurs indépendantes suivant une même loi d’espérance nulle, la moyenne des erreurs divisée par l’écart-type de la loi tend vers une loi normale centrée réduite. De manière plus déterministe, on a l’inégalité de Bienaymé-Tchebyshev

P(|X|>α) ≤ 
nσ2
α2

X est la variable aléatoire somme des n erreurs, α l’erreur et nσ2 la variance de la somme n erreurs supposées indépendantes, cette probabilité tend vers 0 pour n grand si α est d’ordre n, et ne tend pas vers 0 si α est de l’ordre de √n. Exemple : somme de n=400 nombres répartis sur [−1,1] selon la loi uniforme (représentant des erreurs), on divise par √n=20, on effectue plusieurs tirages (par exemple 500) on trace l’histogramme et on compare avec la loi normale de moyenne nulle (l’espérance de la somme) et d’écart-type celui de la loi uniforme.

Il est d’ailleurs souvent trop difficile de calculer une majoration rigoureuse de l’erreur pour des calculs sauf les plus simples. Lorsqu’on doute de la précision d’un calcul, un test peu couteux consiste à refaire ce calcul en utilisant des flottants en précision plus grande et tester si le résultat varie en fonction du nombre de chiffres significatifs utilisés, ou faire varier légèrement les données et observer la sensibilité du résultat. Si on veut travailler en toute rigueur sans pour autant calculer les erreurs à priori, il faut utiliser un logiciel utilisant des intervalles pour représenter les réels (section suivante)

2.3  L’arithmétique d’intervalle.

Certains systèmes de calcul formel peuvent manipuler directement des intervalles réels, par exemple par l’intermédiaire de la bibliothèque C MPFI. Les opérations arithmétiques sur des intervalles renvoient alors le meilleur intervalle possible contenant toutes les valeurs possibles lorsque les opérandes parcourent leurs intervalles respectifs. Exemple en Xcas (version 1.1.1 et ultérieures) : [-1..2]*[-1..2] renvoie [-2..4]. Attention ici on parcourt toutes les valeurs possibles de xy, x ∈ [−1,2], y ∈ [−1,2]. Ce qui est différent du carré d’un intervalle ou plus généralement de l’évaluation d’un polynôme en un intervalle, horner(x^2,[-1..2]) renvoie ainsi [0..4].

Les fonctions disponibles sont souvent moins riches qu’en arithmétique flottante, le calcul d’une fonction non monotone sur un intervalle peut s’avérer délicat, alors que si la fonction est monotone, il suffit de calculer l’image des deux bornes de l’intervalle. Pour les polynômes, Xcas décompose les coefficients en deux parties P=P+P en fonction du signe, puis utilise la monotonie de P+ et P sur ℝ+ et ℝ respectivement.

L’arithmétique d’intervalle dans ℂ est beaucoup plus difficile à mettre en oeuvre puisqu’il n’y a plus d’ordre ni de monotonie, on doit alors s’en remettre à des estimations sur les parties réelles et imaginaires qui ne tiendront pas compte du phénomène ci-dessus sur la différence entre xy, x ∈ [−1,2], y ∈ [−1,2] et x2, x ∈ [−1,2].

2.4  Types composés.

Après les nombres réels, on passe aux nombres complexes : on utilise un couple (partie réelle, imaginaire) de fractions (exacts) ou de flottants et les règles habituelles sur les complexes.

On peut représenter les polynôme de la facon suivante :

Algorithmes de base sur les polynômes : l’évaluation en un point (Horner, cf. TD/TP), la multiplication et division euclidienne et le PGCD (même algorithme que pour les entiers mais avec la division euclidienne des polynômes, attention toutefois aux erreurs d’arrondis)

Les symboles ou noms de variable désignent par exemple le nom d’une inconnue dans un polynôme, ils sont représentés par une chaine de caractére et peuvent être affectés à une valeur pendant une session (la valeur dépend d’un contexte d’exécution et le remplacement du symbole par sa valeur affectée est appelé évaluation).

Les expressions, par exemple sin(x)+2*x^2, elles peuvent être représentées par des arbres. L’évaluation d’une expression consiste à remplacer les symboles de l’expression par leur valeur, puis à effectuer les opérations en tenant compte de la substitution. Il est parfois souhaitable de ne pas effectuer certaines opérations de substitution, on empêche l’évaluation, explicitement ('') ou implicitement (par exemple l’affectation n’évalue pas le symbole qu’on va affecter).

Les fonctions ne doivent pas être confondues avec les expressions, elles associent à leurs arguments une expression. Par exemple sin est une fonction, alors que sin(x) est une expression.

Les conteneurs contiennent plusieurs objets et permettent d’associer à un indice un objet. Il en existe de plusieurs types, par exemple les listes et les séquences dont l’indice est un entier compris entre 1 (ou 0) et la taille (-1), les tables dont l’indice est plus général, et les tableaux (utilisés pour les vecteurs, matrices) qui sont essentiellement des listes ou des listes de listes de même taille. Les séquences sont des listes d’objets ordonnés “non récursifs” (ils ne peuvent contenir des séquences), alors que les listes peuvent contenir des listes, sinon il n’y a pas de différences. Dans les logiciels de calcul formel, la plupart du temps les séquences se notent en indiquant les éléments séparés par des virgules. Les listes s’en distinguent par les délimiteurs []. Il faut prendre garde au fait qu’en général affecter par exemple l[1]:=3; à une variable libre l crée une table et non une liste. Remarque: certains logiciels accédent à certains types de conteneurs uniquement par référence (par exemple maple pour les vecteurs et matrices), dans ce dernier cas une seule copie des objets du conteneur existe si on copie de la manière habituelle une variable contenant un vecteur ou une matrice dans une autre variable, la modification d’un élément du conteneur modifie alors toutes les copies pointant sur ce conteneur. Cette méthode est plus efficace mais peut être surprenante.

3  Algèbre linéaire

3.1  Le pivot de Gauss

Cet algorithme permet de créer des zéros en effectuant des manipulations réversibles sur les lignes d’une matrice. Ces lignes peuvent représenter les coefficients d’un système linéaire, on obtient alors un système linéaire équivalent, ou les coordonnées d’un système de vecteur, on obtient alors les coordonnées d’un système de vecteur engendrant le même sous-espace vectoriel. On peut également représenter 2 matrices A et B reliés par une relation Ax=B, cette relation reste alors vraie au cours et donc après la réduction.

3.1.1  L’algorithme

L’algorithme est le suivant:

  1. on initialise c=1 et l=1, c désigne le numéro de colonne c à réduire, et l le numéro de ligne à partir duquel on cherche un “pivot” (au début l et c valent donc 1, en général les 2 augmentent de 1 à chaque itération)
  2. Si c ou l est plus grand que le nombre de colonnes ou de lignes on arrête.
  3. Si la colonne c n’a que des coefficients nuls à partir de la ligne l, on incrémente le numéro de colonne c de 1 et on passe à l’étape 2. Sinon, on cherche la ligne dont le coefficient est en valeur absolue le plus grand possible (en calcul approché) ou le plus simple possible (en calcul exact), on échange cette ligne avec la ligne l. Puis on effectue pour toutes les lignes sauf l ou pour toutes les lignes à partir de l+1 (selon qu’il s’agit d’une réduction de Gauss complète ou d’une réduction de Gauss sous-diagonale) la manipulation réversible
    Lj ←  Lj −
    ajc
    alc
     Ll 
    On incrémente c et l de 1 et on passe à l’étape 2.

3.1.2  Efficacité de l’algorithme

Si la matrice possède L lignes et C colonnes, le nombre maximal d’opérations pour réduire une ligne est C opérations (une opération=1 multiplication + 1 soustraction, en calculant une seule fois le quotient). Si on a déjà réduit k−1 colonnes, ce nombre est en fait C−(k−1) puisqu’on ne calcule pas les 0 des colonnes déjà réduites. Il y a L−1 lignes à réduire à chaque étape et min(L,C) étapes à effectuer, on en déduit que le nombre maximal d’opérations pour réduire une matrice est (L−1)∑k=1min(L,C)(C−(k−1)). Pour une matrice carrée de taille n, cela fait (n−1) ∑k=1n (n−(k−1)) =n3/2+O(n2) opérations.

Si on effectue une réduction en-dessous de la diagonale, toujours pour une matrice carrée, on trouve ∑k=1n (n−(k−1))(nk)=n3/3+O(n2).

3.1.3  Erreurs d’arrondis du pivot de Gauss

Comme |ajc| ≤ |alc|, une étape de réduction multiplie au plus l’erreur absolue des coefficients par 2. Donc la réduction complète d’une matrice peut multiplier au pire l’erreur absolue sur les coefficients par 2n (où n est le nombre d’étapes de réduction, inférieur au plus petit du nombre de lignes et de colonnes). Ceci signifie qu’avec la précision d’un double, on peut au pire perdre toute précision pour des matrices pas si grandes que ça (n=52). Heureusement, il semble qu’en pratique, l’erreur absolue ne soit que très rarement multipliée par un facteur supérieur à 10, par exemple si on multiplie une matrice n,n à coefficients flottants aléatoires dans [0,1] selon la loi uniforme par son inverse, on trouve une erreur absolue très faible même pour des valeurs de n largement plus grandes que 100.

Par contre, si on ne prend pas la précaution de choisir le pivot de norme maximale dans la colonne, les erreurs d’arrondis se comportent de manière bien moins bonnes, comme le montre l’ :
Exemple
Soit à résoudre le système linéaire

є x + 1.0 y = 1.0 ,    x + 2.0 y = 3.0 

avec є =2−54 (pour une machine utilisant des doubles pour les calculs en flottant, plus généralement on choisira є tel que (1.0+3є)−1.0 soit indistinguable de 0.0).
Si on résoud le système exactement, on obtient x=1/(1−2є) (environ 1) et y=(1−3є)/(1−2є) (environ 1). Supposons que l’on n’utilise pas la stratégie du pivot partiel, on prend alors comme pivot є, donc on effectue la manipulation de ligne L2L2 − 1/є L1 ce qui donne comme 2ème équation (2.0−1.0/є)y=3.0−1.0/є. Comme les calculs sont numériques, et à cause des erreurs d’arrondis, cette 2ème équation sera remplacée par (−1.0/є)y=−1.0/є d’où y=1.0, qui sera remplacé dans la 1ère équation, donnant є x = 1.0−1.0y=0.0 donc x=0.0.
Inversement, si on utilise la stratégie du pivot partiel, alors on doit échanger les 2 équations L2′=L1 et L1′=L2 puis on effectue L2L2′ − є L1′, ce qui donne (1.0−2.0є) y = 1.0 − 3.0 є , remplacée en raison des erreurs d’arrondi par 1.0*y=1.0 donc y=1.0, puis on remplace y dans L1′ ce qui donne x=3.0−2.0y=1.0.
On observe dans les deux cas que la valeur de y est proche de la valeur exacte, mais la valeur de x dans le premier cas est grossièrement eloignée de la valeur correcte.

On peut aussi s’intéresser à la sensibilité de la solution d’un système linéaire à des variations de son second membre. Cela fait intervenir le nombre de conditionnement de la matrice A (voir plus bas) du système (qui est essentiellement la valeur absolue du rapport de la valeur propre la plus grande sur la valeur propre la plus petite), plus ce nombre est grand, plus la solution variera (donc plus on perd en précision).

3.2  Applications de Gauss

3.2.1  Base d’un sous-espace

On réduit la matrice des vecteurs écrits en ligne, puis on prend les lignes non nulles, dont les vecteurs forment une base du sous-espace vectoriel engendré par les vecteurs du départ.

Exemple : base du sous-espac engendré par (1,2,3), (4,5,6), (7,8,9). On réduit la matrice, la 3ème ligne est nulle donc on ne garde que les 2 premières lignes (1,0,−1), (0,1,2) (remarque: une réduction sous-diagonale aurait suffi).

3.2.2  Déterminant

On réduit la matrice (carrée) en notant le nombre d’inversions de ligne, et on fait le produit des coefficients diagonaux, on change le signe si le nombre d’inversions de lignes est impair.

3.2.3  Réduction sous forme échelonnée (rref)

On réduit la matrice puis on divise chaque ligne par son premier coefficient non nul. Si la matrice représentait un système linéaire inversible on obtient la matrice identité sur les colonnes sauf la dernière et la solution en lisant la dernière colonne. La relation conservée est en effet Ax=bx est la solution de l’équation, et à la fin de la réduction A=I.

Par exemple pour résoudre le système





  x + 2y + 3z=
  4x + 5y + 6z=
  7x + 8y=1

on réduit sous forme échelonnée la matrice
ce qui donne [[1,0,0,-9],[0,1,0,8],[0,0,1,-2/3]], d’où la solution x=−9, y=8, z=−2/3.

3.2.4  Inverse

On accolle la matrice identité à droite de la matrice à inverser. On effectue la réduction sous forme échelonnée, on doit obtenir à droite l’identité si la matrice est inversible, on a alors à gauche la matrice inverse. La relation conservée est en effet A x=Bx est l’inverse de la matrice de départ, et en fin de réduction A=I.

Par exemple, pour calculer l’inverse de [[1,2,3],[4,5,6],[7,8,0]], on réduit :

3.2.5  Noyau

On réduit la matrice sous forme échelonnée. Puis on introduit des lignes de 0 afin que les 1 en tête de ligne se trouvent sur la diagonale de la matrice. On enlève ou on rajoute des lignes de 0 à la fin pour obtenir une matrice carrée. Une base du noyau est alors formée en prenant chaque colonne correspondant à un 0 sur la diagonale, en remplaçant ce 0 par -1. On vérifie qu’on obtient bien 0 en faisant le produit de ce vecteur par la matrice réduite. De plus les vecteurs créés sont clairement linéairement indépendants (puisqu’ils sont échelonnés), et il y en a le bon nombre (théorème noyau-image).

Exemple : calcul du noyau de [[1,2,3,4],[1,2,7,12]], on réduit la matrice avec rref, ce qui donne [[1,2,0,-2],[0,0,1,2]], on ajoute une ligne de 0 entre ces 2 lignes pour mettre le 1 de la 2ème ligne sur la diagonale ce qui donne [[1,2,0,-2],[0,0,0,0],[0,0,1,2]], puis on ajoute une ligne de 0 à la fin pour rendre la matrice carrée. On obtient ainsi le système équivalent de matrice :
[[1,2,0,-2],[0,0,0,0],[0,0,1,2],[0,0,0,0]]
La 2ème colonne donne le premier vecteur de la base du noyau, (2,−1,0,0), la 4ème colonne donne le deuxième vecteur (−2,0,2,−1), on vérifie aisément que ces 2 vecteurs forment une famille libre du noyau, donc une base car la dimension du noyau est égale à 2 (dimension de l’espace de départ moins le rang de la matrice, c’est-à-dire le nombre de lignes non nulles de la matrice réduite).

3.3  La méthode de factorisation LU

Dans sa forme la plus simple, elle permet d’écrire une matrice A comme produit de deux matrices triangulaire inférieures et supérieures, ce qui ramène la résolution de système à la résolution de deux systèmes triangulaires. Pour tenir compte d’éléments diagonaux nuls et pour optimiser les erreurs d’arrondi, il est nécessaire d’effectuer des permutations sur les lignes de la matrice.

3.3.1  Interprétation matricielle du pivot de Gauss

On notera l et c le nombre de lignes et colonnes de A (pour éviter la confusion avec le facteur L) et on supposera A non singulière pour simplifier l’exposition.

Lorsqu’on réduit la colonne j d’une matrice à (partiellement réduite) à partir de la ligne j+1 (en supposant Ãj,j ≠ 0), cela revient à multiplier à à gauche par une matrice Lj créée en partant de la matrice identité de taille l où on remplace les 0 colonne j, lignes j+1 à l par le coefficient de la combinaison de ligne effectuée :

li → li − 
Ãi,j
Ãj,j
 lj 

donc :

Lj = 













1...000...0
..................0
0...10......0
0...01......
0...0
− 
Ãj+1,j
Ãj,j
1...
 0...0....0...
0...0
− 
Ãl,j
Ãj,j
0...1
 













On vérifie facilement que l’inverse de cette matrice est

Lj = Lj−1 = 













1...000...0
..................0
0...10......0
0...01......
0...0
Ãj+1,j
Ãj,j
1...
 0...0....0...
0...0
Ãl,j
Ãj,j
0...1
 













Donc A est le produit des matrices Lj par une matrice réduite U qui est triangulaire supérieure

AL1...Ll−1 U 

On vérifie ensuite que le produit des matrices L1...Ll−1 revient à remplacer les coefficients de la colonne j sous la diagonale par ceux de Lj, ce qui donne une matrice L triangulaire inférieure (avec des 1 sur la diagonale). Pour l’obtenir il suffit au cours de l’algorithme de réduction sous-diagonale du pivot de Gauss de stocker le coefficient de la combinaison linéaire dans une matrice initialisée à la matrice identité (on peut aussi le faire en place dans la matrice á réduire).

Attention, le produit Ll−1...L1 ne s’obtient pas en copiant la colonne j de Lj pour j variant de 1 à l−1! On peut l’obtenir en faisant une réduction sous-diagonale de la matrice bloc obtenue en collant A avec la matrice identité ayant l lignes.

3.3.2  Factorisation PA=LU

Si on veut mettre en oeuvre la stratégie du pivot partiel (ou en calcul exact si un coefficient diagonal est nul), il est nécessaire d’intervertir une ligne de la matrice partiellement réduite avec une ligne en-dessous. Cela revient à réduire la matrice A de départ après échange de ces mêmes lignes. En conséquence ce n’est pas A qui est le produit LU mais une matrice obtenue par permutations de lignes de A, que l’on peut écrire comme produit à gauche de A par une matrice de permutation P.

Remarque : si à une étape de réduction, tous les coefficients de la colonne j à partir de la ligne j sont nuls, on peut simplement ignorer cette colonne, on remplit la colonne j de la matrice Lj avec la colonne j+1 de A à partir de la ligne j.

3.3.3  Applications

On peut résoudre des systèmes linéaires par la factorisation LU. En effet soit à résoudre Ax=b. On effectue la permutation de lignes sur A et b (correspondant à la matrice de permutation P), ce qui donne PAx=Pb=LUx, puis on résoud Ly=Pb (système triangulaire inférieur), puis on résoud Ux=y (système triangulaire supérieur).

Comparaison avec la réduction complète sous forme échelonnée de (A|b) :

On peut inverser une matrice en utilisant la décomposition LU. Supposons pour simplifier que la permutation est l’identité. On calcule d’abord L−1 en utilisant le fait que L est triangulaire inférieure, voici comment cela est implémenté dans Xcas (L est noté l) :

first step compute l^-1, 
solve l*a=y for y a canonical basis vector
  a0=y0, a1=y1-l_{1,0}*a0, ..., ak=yk-sum_{j=0..k-1}(l_kj*aj)
if y=(0,..,0,1,0,...0) (1 at position i), 
  a0=..=a_{i-1}=0, a_i=1 and we start at equation k=i+1 and sum_{j=i...}
-> n^3/6 operations
To store the result in place of l
we first compute all the a2 (there is only 1), then all the a3 (2), etc.
a0=y0, a1=y1-l_{1,0}*a0, ..., ak=yk-sum_{j=0..k-1}(l_kj*aj)

Puis on résoud U A−1=L−1 colonne par colonne

second step, solve u*inverse=l^-1 (now under the diagonal)
we compute a column of inverse by solving the system: 
u*col(inverse)=corresponding row of l^-1, 
and overwrite the row of l^-1 by solution
u*[x0,..,xn-1]=[a0,...,an]
x_{n-1}=a_{n-1}/u_{n-1,n-1}
x_{n-2}=(a_{n-2}-u_{n-2,n-1}*x_{n-1})/u_{n-2,n-2}
...
x_k=(a_{k}-sum_{j=k+1..n-1} u_{k,j}x_j)/u_{k,k}
-> n^3/2 operations
To store the solution in place, we first compute all the x_{n-1}
put them in the last line of m, then all the x_{n-2}, etc.

3.4  La factorisation de Cholesky

Dans le cas où la matrice est réelle symétrique définie positive (ou plus généralement hermitienne), on peut obtenir une écriture analogue mais où U est la transconjuguée de L

A=U* U = L L* 

L reste triangulaire inférieure, mais n’a plus des 1 sur sa diagonale en général. Si A est définie positive, on peut rendre l’écriture unique en imposant aux coefficients diagonaux de L d’être réels positifs.

L’algorithme de calcul de U est la traduction matricielle de l’algorithme de Gauss de réduction des formes quadratiques. On a en effet

x* A x=x* U* U x = ||Ux||2 

les lignes de U (ou les colonnes de L) sont donc les coefficients des formes linéaires indépendantes qui interviennent dans l’écriture de la forme quadratique comme somme/différence de carrés de formes linéaires. Si A est définie positive, seules des sommes interviennent, et les variables s’éliminent l’une après l’autre (le coefficient de x2 est forcément non nul lorsqu’on a éliminé x1 et ainsi de suite), ceci explique la forme triangulaire de U et L.

Le calcul de L se fait donc colonne par colonne, en calculant d’abord le coefficient diagonal comme racine carrée du coefficient diagonal αj=√Aj,j. Ensuite on effectue les combinaisons de ligne sous la forme

lj → 
1
αj
lj,     li → αj li − 
Ai,j
αj
 lj 

On peut aussi tout simplement effectuer le produit de L L* et chercher les inconnues en commencant par l1,1 puis on calcule les li,1 pour i>1, etc. En suivant wikipedia, pour une matrice réelle :

L=




l110
l21l220
ln1ln2lnn




aij=
LLT
ij=
n
k=1
likljk
=
min

i,j

k=1
likljk
, 1≤ i,j≤ n 

La matrice A étant symétrique, il suffit que les relations ci-dessus soient vérifiées pour ij, c’est-à-dire que les éléments li,j de la matrice L doivent satisfaire

aij=
i
k=1
likljk
, 1≤ i≤ j≤ n 

Pour i=1, on détermine la première colonne de L

a11=l11l11,    a1j=l11lj1

donc

l11=
a11
,    lj1=
a1j
l11
 (pour  j>1) 

On détermine la i-ième colonne de L (2≤ in) après avoir calculé les i−1 premières colonnes

aii=li1li1+…+liilii,    aij=li1lj1+…+liilji 

d’où

lii
aii
i−1
k=1
lik2
,    lji=
aij
i−1
k=1
likljk
lii
  (pour  j > i)

Pour une matrice hermitienne complexe, il suffit de remplacer likljk par likljk et lik2 par |lik|2.

Le nombre d’opérations arithmétiques à effectuer est asymptotiquement 2 fois plus faible que celui pour LU. En effet, pour la première ligne, il faut 1 racine et n−1 divisions, pour la deuxième ligne, 1 racine, n−1 additions, multiplications et n−2 divisions, ..., pour la i-ième ligne 1 racine, (i−1)(ni) additions, multiplications et n−2 divisions, au final le cout est dominé par les additions et multiplications en n3/6 pour chaque, contre n3/3 pour la factorisation LU. Mais le temps réel dépend de l’optimisation de l’implémentation, en particulier des accès mémoire plus couteux que les opérations arithmétiques, ainsi avec Xcas, l’algorithme de Cholesky n’est pas plus rapide que la factorisation LU.

La commande Xcas correspondante est cholesky et renvoie la matrice L.

3.5  Conditionnement

Le conditionnement mesure la sensibilité de la solution renvoyée d’un système linéaire aux données du problème.

3.5.1  Rappel sur les normes matricielles

On utilisera les normes suivantes sur l’espace vectoriel des matrices n,n :

Si on prend comme norme la norme L2, le calcul de |||A||| nécessite de maximiser √<Ab|Ab> pour b de norme 1, ce qui revient à maximiser √<b|A* A b>. En diagonalisant la matrice hermitienne A* A, on voit qu’il suffit d’en trouver la plus grande valeur propre et d’en prendre la racine carrée. Les valeurs propres de A*A sont appelées valeurs singulières de A (ce sont des réels positifs).

Le calcul de la normes triple infinie est plus simple, si ||x||=1

|(Ax)i| = |
 
j
 ai,j xj| ≤ 
 
j
 |ai,j

avec égalité si les xj valent ± 1 avec le bon choix de signe. On en déduit que

||| A |||=maxi 
 
j
 |ai,j

aussi appelé norme de ligne (on prend le max sur les lignes de la somme des valeurs absolues des coefficients de la ligne).

Pour la norme 1, on trouve la norme colonne (max sur les colonnes de la somme des valeurs absolues des coefficients d’une colonne) :

 
i
 |(Ax)i| = 
 
i
 |
 
j
 ai,j xj| ≤ 
 
j
 |xj
 
i
 |ai,j| ≤ maxj  
 
i
 |ai,j

Pour réaliser la norme triple, on prend xi=0 si ij qui réalise le max, et xj=1.

3.5.2  Nombre de condition

Soit le système linéaire Ax=b de solution x=A−1b, supposons b connu avec une erreur e, alors la solution renvoyée sera x+A−1e, on a donc une erreur relative sur la solution de

||A−1e||
||A−1b||
 = 
||A−1e||
||e||
 
||e||
||b||
||b||
||A−1b||
  ≤ |||A−1|||  
|e|
|b|
 |||A|||  

(la dernière inégalité s’obtient en écrivant b=A(A−1b)). On en déduit que le rapport de l’erreur relative sur la solution par l’erreur relative du second membre est majorée par le produit de la norme de A (en tant qu’application linéaire) par la norme de A−1, ce produit est appelé conditionnement de la matrice A (ou parfois nombre de condition de A en adoptant la terminologie anglo-saxonne).

On remarquera que le conditionnement dépend du choix de la norme sur l’espace vectoriel. Pour la norme 2, on sait que |||A||| est la plus grande valeur singulière de A le même calcul pour A−1 (dont les valeurs singulières sont les inverses des valeurs singulières de A) nous donne alors le :

Théorème 2   Lorsqu’on résoud un système linéaire Ax=b, A matrice connue précisément et inversible, b connu avec une erreur relative en norme L2, l’erreur relative en norme L2 sur x est au plus multipliée par
K2(A)= 
λn
λ1
 
λn [resp. λ1] est la plus grande [resp. plus petite] valeur singulière de A (racines carrées des valeurs propres de A*A).
Ce facteur d’amplification des erreurs relatives est évidemment supérieur ou égal à 1. Il est égal à 1 si la matrice est unitaire (puisque
A est une matrice d’isométrie ou car AA*=I). S’il est de l’ordre de 2c on perdra (au plus) c bits de précision sur la mantisse de x.

Avec Xcas, les valeurs singulières s’obtiennent par l’instruction SVL(A), le conditionnement L2 par cond(A,2) (cond(A,1) donne le conditionnement L1 et cond(A,inf) donne le conditionnement L). Attention, les valeurs singulières de A ne sont pas les valeurs absolues des valeurs propres de A (c’est le cas si A commute avec sa transconjuguée mais ce n’est pas général). On peut utiliser la méthode de la puissance (cf. infra) pour estimer rapidement la plus grande valeur singulière de A (donc sans diagonaliser complètement la matrice A*A), et de même sur A−1 (en utilisant LU ou Cholesky pour trouver les itérées sans calculer A−1).

On remarque que si A est une matrice orthogonale, alors son nombre de conditionnement vaut 1, si on doit utiliser des matrices orthogonales dans un algorithme, ce sera stable numériquement.

On peut aussi prendre la norme L1 sur l’espace vectoriel, dans ce cas la norme de matrice correspondante est la norme de colonne (exercice!), le maximum des sommes valeurs absolues des éléments des colonnes (colNorm(A) en Xcas) et le conditonnement est le produit de colNorm(A) par colNorm(inv(A)) qui est renvoyé par COND(A) en Xcas.

Si la matrice du système A (de nombre de condition noté κ(A)) est elle-même connue avec une certaine incertitude, alors pour ||Δ A|| suffisamment petit, la solution de (AA)(xx)=bb vérifie

|Δ x |
|x|
κ(A)
1−κ(A)
||Δ A||
||A||



|Δ b|
|b|
+
||Δ A||
||A||



En effet, on a

A Δ x =Δ b − Δ A(x+Δ x) ⇒ Δ x = A−1(Δ b − Δ A(x+Δ x))  

donc en norme

|| Δ x || ≤ |||A−1||| (
||Δ b||
||b||
||Ax|| + 
||| Δ A |||
|||A|||
 |||A||| (||x||+||Δ x||))

puis :

|| Δ x ||
||x
 ≤ κ(A) (
||Δ b||
||b||
||| Δ A |||
|||A|||
 (1+
|| Δ x ||
||x

3.6  Quelques méthodes alternatives au pivot

3.6.1  Factorisation QR

La factorisation QR consiste à écrire une matrice A comme produit d’une matrice orthogonale (ou unitaire dans ℂ) et d’une matrice triangulaire supérieure. Les matrices orthogonales ayant un conditionnement de 1 (en norme L2) cette factorisation peut s’obtenir de manière numériquement stable.

Il existe plusieurs algorithmes pour effectuer cette factorisation. On peut voir cette factorisation comme l’orthonormalisation de Gram-Schmidt appliqué aux vecteurs colonnes de A si A est non singulière. Mais la procédure de Gram-Schmidt n’est pas numériquement stable (car on retranche d’un vecteur une combinaison linéaire des vecteurs précdents pour rendre le nouveau vecteur orthogonal et les erreurs d’arrondi s’accumulent rendant l’orthogonalité plus aléatoire). La méthode de Householder utilise des matrices de symétrie par rapport à un hyperplan et ne souffre pas de ce problème d’instabilité. Pour annuler les coefficients de la première colonne c1 de A, on construit le vecteur u=c1 ± ||c1|| e1e1 est le premier vecteur de base et le signe ± est le signe de la première composante de c1 (pour assurer la stabilité numérique). On fait alors la symétrie Q1 par rapport à l’hyperplan H orthogonal à u, qui laisse H invariant et transforme u en −u. Comme u est vecteur directeur de la bissectrice intérieure ou extérieure de c1 et e1, la symétrie échange ces deux vecteurs, éventuellement au signe près. La matrice A1=Q1A a donc comme première colonne un multiple de e1, on continue ensuite en faisant le même raisonnement sur la matrice A1 en se limitant et lignes et colonnes d’indice ≥ 2. Après n−1 itérations, on a An−1=Qn−1...Q1A qui est triangulaire supérieure d’où la factorisation annoncée.

Matriciellement, Q1=I−2vv*v=u/||u||, pour calculer Q1A il faut effectuer A−2vv*A, on calcule donc w=v*A en n2 opérations (ou une opération est une addition et une multiplication) puis on soustrait (2v)i wj de aij en n2 opérations. En faisant de même aux étapes qui suivent, sans tenir compte de la simplification progressive du vecteur v, on effectue 2n3 opérations. La constante 2 peut être un peu améliorée en tenant compte des 0 initiaux de v aux étapes 2 et ultérieures, elle est toutefois supérieure à LU (et Cholesky), mais en contrepartie la méthode est très stable numériquement.

On peut aussi utiliser des rotations (méthode de Givens) pour annuler les coefficients de A sous la diagonale. C’est par exemple très efficace pour des matrices tridiagonales.

Lorsque la matrice A n’est pas carrée, mais possède n lignes et c colonnes, la factorisation QR est encore possible, la matrice Q est alors une matrice carrée d’ordre n, et R a les mêmes dimensions que A. En particulier si A a plus de lignes que de colonnes (matrice verticale), R aussi. Par exemple si A est de rang maximal c, alors R se décompose en un premier bloc c,c inversible R1 et un deuxième bloc nc,c entièrement nul.

Applications :
On peut alors écrire Ax=b sous la forme QRx=b donc Rx=Q*b qui est un système triangulaire supérieur, donc résoudre Ax=b en O(n2) opérations une fois la factorisation QR effectuée.

Mais c’est surtout pour résoudre au sens des moindres carrés un système sur-déterminé que la factorisation QR trouve tout son intérêt. Soit A la matrice d’un système sur-déterminé avec n lignes et c colonnes, n>c (matrice “verticale” ou mince). Le système n’a en général pas de solution, on cherche alors à minimiser || Axb ||22. Ceci revient à chercher la projection orthogonale de b sur Im(A). Soit Ax cette projection, on a b=Ax+p avec p orthogonal à Im(A) donc dans Ker(A*), donc A*b=A*(Ax+p)=A*Ax.

Proposition 3   La solution x du problème de minimisation de || Axb ||2 est donnée par A*Ax=A*b

Si A=QR on a alors R* R x=R* Q*b. Si A est de rang maximal c, on décompose R en (R2R1) avec R1 inversible de taille c et R2=0, donc R*=(R1*, 0) et le système devient R1*R1x =(R1*,0) Q*b et comme R1 est inversible, on a R1x= les c premières lignes de Q*b. Il vaut mieux résoudre ce système que A* Ax=A*b car le conditionnement de A* A est le carré du conditionnement de A (le nombre d’opérations est un peu supérieur : QR a une constante nettement plus grande que Cholesky mais il n’y a pas de multiplication de matrice à faire, et le résultat est plus précis).

Notons qu’on peut résoudre de manière assez analogue un système sous-déterminé au sens des moindres carrés. On considère toujours le système Ax=b mais cette fois-ci A est une matrice “horizontale” (nombre de lignes n<c le nombre de colonnes). Si b ∈ Im(A) (on peut génériquement supposer que A est de rang n), il y a une infinité de solutions, on cherche alors la solution de norme minimale. On cherche donc le projeté orthogonal de l’espace affine des solutions (dirigé selon Ker(A)) sur l’orthogonal de Ker(A) qui est ImA*. Ainsi x=A*y vérifie Ax=b, donc AA*y=b. Si A est de rang n alors AA* est inversible et y=(AA*)−1b donc x=A (AA*)−1b. En utilisant la factorisation QR de A* (attention pas celle de A), on a AA*=R*R, donc x=QR(R*R)−1b. Si A est de rang n, on a alors en posant R=(0R1) avec R1 matrice n,n inversible :

x=Q

R1
0


R1−1 R1*−1b=Q

R1*−1 
0


b=Q1 R1*−1 b

Q1 désigne les n premières colonnes de Q.

Proposition 4   Soit A une matrice ayant n lignes et c>n colonnes. Si A est de rang maximal n, posons A*=QR, R1 les n premières lignes de R1, Q1 les n premières colonnes de Q. La solution de Ax=b de norme minimale est donnée par Q1cc est la solution de R1*c=b (système triangulaire).

Exemple 1 : la régression linéaire :
On se donne n points de coordonnées (xi,yi) où les xi sont distincts 2 à 2 et on cherche la droite qui approche le mieux ces points au sens de l’écart en y au carré. Il s’agit donc de trouver α et β qui minimisent :

n
j=1
 (α+β xiyi)2 

Ici on a :

A=


1x1 
...... 
1xn  
 


,     x=

α 
β 


,    b=


y1  
...  
yn  
 


Pour trouver α, β de manière précise, on préferera donc calculer la factorisation QR de A plutôt que calculer A*A, R1 est une matrice carrée d’ordre 2, Q une matrice carrée d’ordre n, la première colonne de Q est le vecteur de coordonnées ± (1,...,1)/√n. Exemple : on se donne les points (1,2), (2,1), (3,5), (4,4). La matrice A et b sont donc




Exemple 2: régression polynomiale :
C’est une généralisation de l’exemple précédent, on cherche un polynôme P(x) de degré fixé tel que ∑i (yiP(xi))2 soit minimal. Si on choisit la base canonique des polynômes, la matrice A est alors une matrice de Vandermonde et son conditionnement est en général très mauvais. De plus, plus le degré augmente, plus l’extrapolation est hasardeuse :




On préfèrera utiliser une base de polynômes orthogonaux pour le produit scalaire ∑i P(xi)Q(xi) ou minimiser une autre fonctionnelle quadratique, par exemple une intégrale du carré de la différence entre une fonction à approcher et un polynôme ou un polynôme trigonométrique.

Remarque : La factorisation QR peut s’obtenir en calculant la factorisation de Cholesky de A* A qui donne R, mais on lui préfère une méthode de calcul direct pour des raisons de stabilité, en tout cas pour des matrices denses (pour des matrices creuses des considérations d’efficacité peuvent faire préférer la factorisation de Cholesky).

3.6.2  Jacobi, Gauss-Seidel, relaxation

Lorsqu’on a une matrice creuse (peu d’éléments non nuls), l’algorithme du pivot de Gauss a tendance à densifier rapidement la matrice réduite (surtout avec le pivot partiel où on ne controle pas le nombre de zéros de la ligne contenant le pivot). Il peut alors être intéressant d’utiliser des méthodes alternatives ne faisant intervenir que des produits de matrice, donnant éventuellement un résultat seulement approché.

Par exemple pour calculer l’inverse d’une matrice A=MN, avec M facile à inverser (par exemple diagonale) et N petit en norme et creuse, on peut écrire :

A−1=(MN)−1 = (M(IM−1N))−1=(I+M−1N+(M−1N)2+...)M−1

De même pour résoudre un système linéaire Ax=b avec A=MN, on considère la suite Mxn+1Nxn=b, donc xn+1 est obtenu en résolvant le système :

Mxn+1=b+N xn,    x0=0

pour laquelle on vérifiera les hypothèses du théorème du point fixe, il suffit par exemple de vérifier que la plus grande valeur singulière de M−1N est strictement plus petite que 1. Lorsque la matrice N n’est pas creuse, le procédé est intéressant pour résoudre approximativement un système si n est grand et si chaque itération est en O(n2) (ceci veut dire qu’on ne calcule pas M−1 sauf si c’est évident, par exemple si M est diagonale), mais le procédé n’est pas intéressant pour le calcul de l’inverse de A.

Notons D la partie diagonale de A, L sa partie triangulaire inférieure stricte, U sa partie triangulaire supérieure stricte, La méthode de Jacobi utilise pour M la diagonale D de A, alors que la méthode de Gauss-Seidel prend pour M la partie triangulaire inférieure D+L de A (diagonale comprise). Pour Jacobi, on a donc

xn+1=D−1(b+(DA)xn)=xn+D−1(bAxn)

En Xcas, l’instruction préprogrammée est jacobi_linsolve, on peut aussi programmer la méthode par

jacobi(A,b,eps,N):={
  local D,x0,x1,n,j;
  n:=size(A);
  D:=diag(A).^-1;
  x0:=seq(0,n);
  pour j de 1 jusque N faire
    x1:=b-A*x0;
    si l2norm(x1)<eps alors return x0; fsi;
    x0:=x0+D .*x1;
  fpour;
  return "non convergent";
}:;


Un cas simple où on a convergence :

Proposition 5   Lorsque la matrice A est à diagonale strictement dominante, c’est-à-dire que l’élément diagonal est en valeur absolue strictement supérieur à la somme des éléments non diagonaux de la même ligne :
|aii| > 
 
j≠ i
 |aij
la méthode de Jacobi converge.

En effet, ||M−1Nx||< || x|| car :

|(Nx)i| ≤ 
 
j ≠ i
 |aij| |xj| ≤ 
 
j ≠ i
|aij| ||x||< |aii| || x||    ⇒    |M−1(Nx)i| <  || x||

On retrouve ce cas pour une petite perturbation d’une matrice diagonale, par exemple 








Pour n assez grand, la méthode de Jacobi devient plus rapide. Cela se vérifie encore plus vite si A est une matrice creuse.

Pour Gauss-Seidel, le calcul de M−1 n’est pas effectué, on résoud directement le système triangulaire Mxn+1=b+Nxn soit

(D+L)xn+1=bUxn

Gauss-Seidel est moins adapté à la parallélisation que Jacobi. On adapte le programme précédent

seidel(A,b,N,eps):={
  local L,U,x0,x1,n,j;
  n:=size(A);
  L:=diag(A,left);
  U:=A-L;
  x0:=seq(0.0,n);
  pour j de 1 jusque N faire
    x1:=b-U*x0;
    x1:=linsolve(L,x1);
    si l2norm(x1-x0)<eps*l2norm(x0) alors return x0; fsi;
    x0:=x1;
  fpour;
  return "non convergent";
}:;


Dans la méthode de relaxation, on pose pour M la matrice triangulaire inférieure M=1/ωD+L où ω >0, donc N=(1/ω−1) DU et on utilise la récurrence Mxn+1=b+Nxn donc Mxn+1=b+(MA)xn puis M(xn+1xn)=bAxn puis

(D+ω L) (xn+1xn) = ω(bAxn)

On remarque que Gauss-Seidel correspond à ω=1. L’instruction correspondante de Xcas est gauss_seidel_linsolve et peut prendre en paramètre le coefficient de relaxation, dont la valeur par défaut est 1.

Proposition 6   (Convergence ): si A=MN est une matrice symétrique définie positive et si M*+N est définie positive, alors la méthode converge.

On utilise la norme correspondant à la forme quadratique de matrice A et on calcule la norme subordonnée de M−1N, on a M−1Nx=xy avec y=M−1Ax donc

|| M−1Nx ||A2=<xy|A (xy))> 
 =<x|Ax>+<y|Ay>−<y|Ax>−<x|Ay
 =<x|Ax>+<y|Ay>−<y|My> −<My|y>
 =<x|Ax>−<y|(M*+N)y>
 <<x|Ax>=||x||A2

Conséquence : si A symétrique définie positive, alors Gauss-Seidel converge, car M*+N=D. Pour la relaxation, on a M*+N=(2/ω −1)D qui est définie positive si ω<2.

Remarque : Jacobi et Gauss-Seidel sont implémentées dans les commandes Xcas jacobi_linsolve et gauss_seidel_linsolve.

3.6.3  Le gradient conjugué

Il s’agit de résoudre Ax=b, où A est définie positive. Si on a une base orthogonale pour le produit scalaire induit par A, on peut calculer la j-ième coordonnée de x dans cette base en faisant le produit scalaire de Ax=b par le j-ième vecteur de la base. On construit donc petit à petit une base orthogonale pour A par un procédé à la Gram-Schmidt, mais on ne part pas de la base canonique : on construit cette famille orthogonale pour A en meme temps qu’on calcule les composantes de x. On pose x0=0, à la i-ième itération si Axib=0 on a terminé, sinon ri=Axib est linéairement indépendant des éléments de la famille orthogonale dejà construite, on complète la famille orthogonale avec un nouveau vecteur, on calcule la i+1-ième composante de x sur la famille orthogonale, et on ajoute le tout à xi pour obtenir xi+1. On s’arrête en au plus la dimension itérations lorsque la famille orthogonale est devenue une base.

La commande conjugate_gradient(A,b) de Xcas permet de faire ce calcul, on peut préciser une valeur initiale de recherche x0 et une précision eps en tapant conjugate_gradient(A,b,x0,eps). Voir aussi le menu Exemple, analyse, gradconj

4  Approximation polynomiale

On présente dans cette section quelques méthodes d’approximation de fonctions par des polynômes sur un intervalle, la section suivante présente des méthodes d’approximation près d’un point ou de l’infini.

4.1  Polynôme de Lagrange

Étant donné la facilité de manipulation qu’apportent les polynomes, on peut chercher à approcher une fonction par un polynôme. La méthode la plus naturelle consiste à chercher un polynôme de degré le plus petit possible égal à la fonction en certains points x0,...,xn et à trouver une majoration de la différence entre la fonction et le polynôme. Le polynome interpolateur de Lagrange répond à cette question.

4.1.1  Existence et unicité

Soit donc x0,...,xn des réels distincts et y0,...,yn les valeurs de la fonction à approcher en ces points (on posera yj=f(xj) pour approcher la fonction f). On cherche donc P tel que P(xj)=yi pour j ∈ [0,n].

Commencons par voir s’il y a beaucoup de solutions. Soit P et Q deux solutions distinctes du problème, alors PQ est non nul et va s’annuler en x0, ...,xn donc possède n+1 racines donc est de degré n+1 au moins. Réciproquement, si on ajoute à P un multiple du polynome A=∏j=0n (Xxj), on obtient une autre solution. Toutes les solutions se déduisent donc d’une solution particulière en y ajoutant un polynome de degré au moins n+1 multiple de A.

Nous allons maintenant construire une solution particulière de degré au plus n. Si n=0, on prend P=x0 constant. On procède ensuite par récurrence. Pour construire le polynôme correspondant à x0,...,xn+1 on part du polynoôme Pn correspondant à x0,...,xn et on lui ajoute un multiple réel de A

Pn+1=Pn+ αn+1 
n
j=0
 (Xxj

Ainsi on a toujours Pn+1(xj)=yj pour j=0,..n, on calcule maintenant αn+1 pour que Pn+1(xn+1)=yn+1. En remplacant avec l’expression de Pn+1 ci-dessus, on obtient

Pn(xn+1)+  αn+1 
n
j=0
 (xn+1xj) = yn+1 

Comme tous les xj sont distincts, il existe une solution unique :

αn+1=
yn+1Pn(xn+1)
n
j=0
 (xn+1xj)

On a donc prouvé le :

Théorème 7   Soit n+1 réels distincts x0,...,xn et n+1 réels quelconques y0,...,yn. Il existe un unique polynôme P de degré inférieur ou égal à n, appelé polynome de Lagrange, tel que :
P(xi)=yi

Exemple : déterminons le polynome de degré inférieur ou égal à 2 tel que P(0)=1, P(1)=2, P(2)=1. On commence par P0=1. Puis on pose P1=P0+ α1X=1+ α1X. Comme P(1)=2=1+ α1 on en tire α1=1 donc P1=1+X. Puis on pose P2=P1+ α2X(X−1), on a P2(2)=3+2 α2=1 donc α2=−1, finalement P2=1+XX(X−1).

4.1.2  Majoration de l’erreur d’interpolation.

Reste à estimer l’écart entre une fonction et son polynome interpolateur, on a le :

Théorème 8   Soit f une fonction n+1 fois dérivable sur un intervalle I=[a,b] de , x0,...,xn des réels distincts de I. Soit P le polynome de Lagrange donné par les xj et yj=f(xj). Pour tout réel xI, il existe un réel ξx ∈ [a,b] (qui dépend de x) tel que :
  f(x)−P(x) = 
f[n+1]x)
(n+1)!
 
n
j=0
(xxj)      (1)

Ainsi l’erreur commise dépend d’une majoration de la taille de la dérivée n+1-ième sur l’intervalle, mais aussi de la disposition des points xj par rapport à x. Par exemple si les points xj sont équidistribués, le terme |∏j=0n(xxj)| sera plus grand près du bord de I qu’au centre de I.

Preuve du théorème : Si x est l’un des xj l’égalité est vraie. Soit

C=(f(x)−P(x))/
n
j=0
(xxj

on considère maintenant la fonction :

g(t)=f(t)−P(t) − C 
n
j=0
(txj

elle s’annule en xj pour j variant de 0 à n ainsi qu’en x suite au choix de la constante C, donc g s’annule au moins n+2 fois sur l’intervalle contenant les xj et x, donc g′ s’annule au moins n+1 fois sur ce même intervalle, donc g′′ s’annule au moins n fois, etc. et finalement g[n+1] s’annule une fois au moins sur cet intervalle. Or

g[n+1] = f[n+1] − C (n+1)!

car P est de degré inférieur ou égal à n et ∏j=0n(xxj) − xn+1 est de degré inférieur ou égal à n. Donc il existe bien un réel ξx dans l’intervalle contenant les xj et x tel que

C=
f[n+1]x)
(n+1)!
 







Attention, l’erreur d’interpolation peut devenir très grande lorsqu’on utilise beaucoup de points d’interpolation.

4.1.3  Calcul efficace du polynôme de Lagrange.

Avec la méthode de calcul précédent, on remarque que le polynôme de Lagrange peut s’écrire à la Horner sous la forme :

 P(x)=α0 + α1 (xx0) + ... + αn (xx0)...(xxn−1
 =α0 + (xx0)( α1 + (xx1)(α2 +  ... + (xxn−2)(αn−1+(xxn−1) αn)...))

ce qui permet de le calculer rapidement une fois les αi connus. On observe que

α0=f(x0),    α1=
f(x1)−f(x0)
x1x0
 

On va voir que les αk peuvent aussi se mettre sous forme d’une différence. On définit les différences divisées d’ordre n par récurrence

f[xi]=f(xi),    f[xi,...,xk+i+1]=
f[xi+1,...,xk+i+1]−f[xi,...,xk+i]
xk+i+1xi
 

On va montrer que αk=f[x0,...,xk]. C’est vrai au rang 0, il suffit donc de le montrer au rang k+1 en l’admettant au rang k. Pour cela on observe qu’on peut construire le polynôme d’interpolation en x0,...,xk+1 à partir des polynômes d’interpolation Pk en x0,...,xk et Qk en x1,...,xk+1 par la formule :

Pk+1(x)= 
(xk+1x)Pk + (xx0)Qk
xk+1x0

en effet on vérifie que Pk+1(xi)=f(xi) pour i∈ [1,k] car Pk(xi)=f(xi)=Qk(xi), et pour i=0 et i=k+1, on a aussi Pk+1(x0)=f(x0) et Pk+1(xk+1)=f(xk+1). Or αk+1 est le coefficient dominant de Pk+1 donc c’est la différence du coefficient dominant de Qk et de Pk divisée par xk+1x0, c’est-à-dire la définition de f[x0,...,xk+1] en fonction de f[x1,...,xk+1] et f[x0,...,xk].

Exemple : on reprend P(0)=1, P(1)=2, P(2)=1. On a

xif[xi]f[xi,xi+1]f[x0,x1,x2
0
1
  
  
(2−1)/(1−0)=
1
 
12 
(−1−1)/(2−0)=
-1
   
  (1−2)/(2−1)=−1 
21  

donc P(x)=1+(x−0)(1+(x−1)(−1))=1+x(2−x).

On peut naturellement utiliser l’ordre que l’on souhaite pour les xi, en observant que le coefficient dominant de P ne dépend pas de cet ordre, on en déduit que f[x0,...,xk] est indépendant de l’ordre des xi, on peut donc à partir du tableau ci-dessus écrire P par exemple avec l’ordre 2,1,0, sous la forme

P(x)=1+(x−2)(−1+(x−1)(−1))=1+(x−2)(−x)

Le nombre d’opérations nécessaires pour faire ce calcul est proportionnel à n2.

La commande Xcas interp ou son synonyme lagrange effectue ce calcul. Pour avoir les différences divisées, on peut créer le programme suivant :

dd(X,Y):={ // Algorithme des différences divisées
  local k,l,n,A,old,cur;
  si size(X)!=size(Y) alors return "erreur" fsi;
  n:=size(X)-1;
  A:=[Y[0]];
  old:=Y;
  pour k de 1 jusque n faire
    // calcul de cur en fonction de old
    cur:=[];
    pour l de 0 jusque n-k faire
      cur[l]:=(old[l+1]-old[l])/(X[l+k]-X[l])
    fpour;
    A[k]:=cur[0];
    old:=cur;
  fpour;
  retourne A;
}:;



(N.B. pour rendre ce programme optimal, il faudrait utiliser l’affectation en place =< au lieu de :=)

4.1.4  Sensibilité aux erreurs sur les données.

Si les yj sont connus avec une certaine erreur, alors le polynôme d’interpolation est connu de manière approchée. Plus précisément, si on note

πj(x)=
 
k ≠ j
 
xxk
xjxk

le j-ième polynome de Lagrange valant 1 en xj et 0 ailleurs, l’erreur vaut :

 
j
 (ỹjyj) πj(x

Si l’erreur relative sur les yj est majorée par є, l’erreur sur le polynôme d’interpolation est majorée par :

є   maxj|yj
 
j
 |πj(x)| 

il y a amplification de l’erreur par un facteur majoré par 

maxx ∈ [a,b] 
n
j=0
 |πj(x)| 

Ce facteur s’appelle constante de Lebesgue relative à la subdivision x0,...,xn de [a,b]. On peut le calculer numériquement pour une subdivision équidistribuée, et montrer qu’il croit comme 2n+1/e n ln(n), par exemple pour n=40, il vaut environ 5e9. Illustration avec Xcas :




puis essayer avec n=20. Pour n=40, en observant que le max est atteint dans [0,1], on peut remplacer les valeurs absolues par la bonne puissance de -1

on a alors un polynôme, dont on calcule l’abscisse du maximum par

puis
qui donne environ 4.7e9.

4.2  Interpolation aux points de Tchebyshev

L’idée la plus naturelle pour interpoler un polynôme en n+1 points d’un intervalle [a,b] consiste à couper en n morceaux de même longueur. Mais ce n’est pas le plus efficace car le terme |∏j=0n(xxj)| est plus grand près des bords. Il est donc plus judicieux d’avoir plus de points près des bords et moins à l’intérieur. C’est là qu’interviennent les polynômes de Tchebyshev, ils sont définis par développement de cos(nx) en puissances de cos(x) :

Tn(cos(x))=cos(nx

Sur [−1,1], le polynôme Tn vaut en valeur absolue au plus 1, et atteint cette valeur exactement n+1 fois lorsque x=kπ/n donc X=cos(x)=cos(kπ/n).


De plus cette majoration est optimale. En effet soit U un polynôme de degré au plus n qui vérifie |U|<1 sur [−1,1] et tel que U ait le même coefficient dominant que Tn. Alors la différence TnU est du signe de Tn en X=cos(kπ/n), k ∈ [0,n] puisqu’en ces points Tn est extrêmal de valeur absolue 1. Donc TnU s’annule n fois sur [−1,1], mais son degré est au plus n−1 donc Tn=U absurde.

On a donc intérêt à prendre les abscisses des points d’interpolation en les racines tn de Tn

a+b
2
 + 
ab
2
 tn,    tn=cos((k+
1
2
)
π
n
),     k=0..n−1

On pourra observer que le phénomène de Runge qui apparait par exemple pour f(x)=1/(25x2+1) sur [−1,1] avec des points d’interpolation équidistants n’apparait plus si on prend des points de Tchebyshev.







Ceci est relié à la constante de Lebesgue qui pour des points de Tchebyshev vaut un peu moins de 4 pour n<100 (se comporte comme 2/πln(n) pour n grand), on peut montrer que les polynômes de Lagrange aux points de Tchebyshev convergent uniformément vers 1/(25x2+1) (c’est plus généralement vrai pour toute fonction C1 sur l’intervalle).

Remarque : ce n’est pas le polynôme de meilleure approximation, de f (celui qui minimise la norme L de la différence) car la dérivée n+1-ième varie en général sur [a,b]. Mais il est trop difficile de le calculer en général.

Exemple de calcul explicite de constante de Lebesgue pour n=40 avec Xcas

purge(x);
t(k,n):={ 
 local T; 
 T:=seq(cos(pi*(k+.5)/(n+1)),k,0,n);
 return product((x-T[j])/(T[k]-T[j]),j,0,k-1)*
  product((x-T[j])/(T[k]-T[j]),j,k+1,n);
}:;



4.3  Interpolation de Hermite

Si on fait tendre un des points d’interpolation vers un autre, la donnée de la valeur en ces 2 points serait redondante, elle est remplacée par la valeur de la dérivée. Dans le calcul des différences divisées ci-dessus on fera comme si les 2 points étaient distincts et successifs, disons xi et xi+1, on remplace le rapport indéterminé

f(xi+1)−f(xi)
xi+1xi
 = 
0
0

par f′(xi). On montre qu’une fois ce changement réalisé tout le reste est identique (y compris la majoration d’erreur). On peut bien sur généraliser au cas de plusieurs paires de points identiques ou des multiplicités plus grandes faisant intervenir des dérivées d’ordre supérieures, dans ce cas la différence divisée f[xi,...,xi+m] sera remplacée par f[m](xi)/m!.

4.4  Polynômes de Bernstein et courbes de Bézier

Les polynômes de Bernstein de degré m sont les

Bkn (x) = 

n 
k 
 

xk (1−x)nk

On reconnait la probabilité d’avoir k succès si on effectue n tirages indépendants (avec remise) avec probabilité x ∈ [0,1] de succès par tirage. Ceci donne une relation de récurrence

Bk+1n+1=(1−x)Bk+1n+x Bkn 

qui peut servir à calculer les Bim. On en déduit aussi que l’espérance de k selon cette loi vaut nx (somme de n variables d’espérance x) et l’espérance de (knx)2 vaut nx(1−x) (variance de la somme de n variables indépendantes de variance x). On en déduit qu’on peut approcher uniformément une fonction continue sur un intervalle [a,b] par des polynômes, en se ramenant à a=0, b=1, on pose :

Pn(x) = 
n
k=0
 f(
k
n
Bkn(x

En effet, par continuité uniforme de f sur [0,1], pour є>0, il existe δ>0 tel que |xy|<δ ⇒ |f(x)−f(y)|<є/2, dans

Pn(x)−f(x)=
n
k=0
 (f(
k
n
)−f(x)) Bkn(x)

on décompose la somme sur k en deux parties, |k/nx|<δ et |k/nx| ≥ δ, pour la première somme, on majore |f(k/n)−f(x)| par є/2 puis par ∑k=0n, pour la deuxième somme, on majore par 2|f| et on utilise 1<(k/nx)22=1/n22 (knx)2 pour se ramener au calcul de la variance de k, au final

|Pn(x)−f(x)| ≤ 
є
2
 + 
1
n2 δ2
 nx(1−x)|f|

il suffit de choisir n assez grand pour rendre le membre de droite plus petit que є.

Les polynômes de Bernstein ne sont pas des polynômes interpolateurs aux points k/n,   0<k<n, et la convergence n’est pas forcément très rapide. On les utilise pour approcher rapidement des morceaux de courbes, si on se donne des “points de controle” A0, ..., An on construit la courbe paramétrée

A(t) = 
n
k=0
 Ak 

n
i 
 

xi (1−x)ni 

appelée courbe de Bézier. En pratique on les utilise pour n=3.

4.5  Polynômes orthogonaux.

Il s’agit d’une autre méthode d’approximation, particulièrement important pour l’intégration : les polynômes de meilleur approximation au sens de normes L2 ou L2 à poids w(x)>0 sur l’intervalle de bornes α et β (finis ou infinis). On considère le produit scalaire

<f|g>=
β


α
w(xf(x)g(x)  dx

et on projette alors la fonction à approcher sur une base de polynômes orthogonaux de degrés croissants construit par la procédure de Gram-Schmidt à partir de la base canonique pour le produit scalaire ci-dessus.

Proposition 9   Le polynôme Pn de degré n obtenu par orthogonalisation de Gram-Schmidt pour le produit scalaire αβw(x) f(x)g(x) dx possède n racines réelles.

En effet, soit r le nombre de racines réelles de Pn, on pose f=Pn et g=∏i=1r (xxi), wfg est de signe constant et non identiquement nul donc <f|g> ≠ 0 donc r=n sinon g serait de degré strictement plus petit que n donc orthogonal à Pn.

On peut aussi construire ces polynômes en cherchant les valeurs propres de

T(f)=
1
w
 (awf′)′

a>0 est un polynôme de degré au plus 2 sur ]α,β[, tel que aw s’annule (ou tend vers 0) aux bornes de l’intervalle d’intégration : si α et β sont finis a(x)=(x−α)(β−x). On a alors

<T(f)|g>=
β


α
(awf′)′g = [awfg]αβ
β


α
awgf′ =−
β


α
awgf

car le terme tout intégré s’annule (puisque aw s’annule en α et β). Donc <T(f)|g>=<f|T(g)> est symétrique, les vecteurs propres de T correspondant à des valeurs propres distinctes sont donc orthogonaux entre eux. Pour trouver ces valeurs propres/polynômes vecteurs propres, on écrit T(Pn)=λn Pn pour un polynôme Pn de degré n. Si aw′/w est un polynôme de degré au plus 1, le terme de degré dominant de cette équation donne la valeur de λn et les termes de degré plus petits permettent en général de déterminer de manière unique les coefficients de Pn en fonction du coefficient dominant.

Pour certains poids w(x) standards, les polynômes obtenus ont un nom :

Ainsi, les polynômes de Legendre vérifient

((1−x2)Pn′)′=λn Pn

Le terme de degré n de cette équation donne λn=−n(n+1), le terme sous-dominant est nul. Plus générallement, le terme de degré k vérifie

  (k+2)(k+1)pk+2k(k+1)pk=−n(n+1)pk ⇒ pk+2=
k(k+1)−n(n+1)
(k+2)(k+1)
 pk     (2)

Ceci permet de calculer le polynôme Pn, on normalise Pn par Pn(1)=1.




Le calcul de ces polynômes peut aussi se faire par une récurrence à 2 crans du type :

  Pn+1=(anxbn)PncnPn−1     (3)

La valeur de an est définie par la normalisation de la famille de polynômes, ensuite bn et cn sont déterminés respectivement par

<Pn+1|Pn>=0 ⇒ bn<Pn|Pn>=an<xPn|Pn>

et

<Pn+1|Pn−1>=0 ⇒ cn<Pn−1|Pn−1>=an<xPn|Pn−1>=an<Pn|xPn−1>

donc

cn<Pn−1|Pn−1>=
an
an−1
<Pn|Pn+bn−1Pn−1+cn−1Pn−2> =
an
an−1
<Pn|Pn>

Les autres relations d’orthogonalité <Pn+1|Pj>=0, j<n−1 sont automatiquement vérifiées puisque <xPn|Pj>=<Pn|xPj>=0 et Pn est orthogonal aux polynômes de degré ≤ n−1.

Dans l’exemple des polynômes de Legendre, les polynômes obtenus sont pairs si n est pair et impairs sinon, la relation de récurrence a donc un coefficient bn nul. La convention de normalisation usuelle est Pn(1)=1, on peut montrer qu’on a an=(2n+1)/(n+1) et <Pn|Pn>=2/(2n+1) donc cn=n/(n+1).

(n+1)Pn+1=(2n+1)xPnnPn−1

En effet la normalisation en 1 donne 1=ancn donc cn=an−1, les termes de degré n+1 et n−1 de Pn+1=anxPncnPn−1 donnent

pn+1,n+1=anpn,n,    pn+1,n−1=an pn,n−2−(an−1)pn−1,n−1

On applique alors (2) et on déduit de la deuxième équation ci-dessus :

n(n+1)
2(2n+1)
 pn+1,n+1= −
n(n−1)
2(2n+1)
 an pn,n−(an−1)pn−1,n−1

puis on applique la première :

n(n+1)
2(2n+1)
 −
n(n−1)
2(2n+1)
=
an−1
anan−1

d’où

an=
1
1−
n2
4n2−1
 an−1







Remarques

4.6  Les splines

Il s’agit de fonctions définies par des polynomes de degré borné sur des intervalles, dont on fixe la valeur aux extrémités des intervalles (comme pour le polynome de Lagrange) ce qui rend la fonction continue, de plus on exige un degré de régularité plus grand, par exemple etre de classe C2. Enfin, on fixe des conditions aux bornes de la réunion des intervalles, par exemple avoir certaines dérivées nulles.

Par exemple supposons qu’on se donne n intervalles, donc n+1 points x0,...,xn, on se fixe une régularité Cd−1. Ceci entraine (n−1)d conditions de recollement, on y ajoute n+1 conditions de valeur en x0,...,xn, on a donc nd+1 conditions, la borne sur le degré des polynomes doit donc etre d (ou plus, mais d suffit) ce qui donne n(d+1) degrés de liberté, on peut donc ajouter d−1 conditions, par exemple pour les splines naturelles, on impose que les dérivées d’ordre d/2 à d−1 soient nulles en x0 et xn (si d est pair, on commence à la dérivée d/2+1-ième nulle en xn).

Pour trouver les polynomes, on doit donc résoudre un grand système linéaire. Une méthode permettant de diminuer la taille du système linéaire à résoudre dans le cas des splines naturelles consiste à se fixer n inconnues z0,..,zn−1 représentant les dérivées d-ième de la spline f en x0 sur [x0,x1] à xn−1 sur [xn−1,xn], et (d−1)/2 inconnues fj, représentant la valeur de la dérivée de f en x0 pour j variant de 1 à (d−1)/2. On peut alors écrire le polynome sur l’intervalle [x0,x1] car on connait son développement de Taylor en x0. On effectue un changement d’origine (par application répétée de Horner) en x1. On obtient alors le polynome sur [x1,x2] en remplaçant uniquement la dérivée d-ième par z1. On continue ainsi jusqu’en xn−1. Le système s’obtient en calculant la valeur du polynome en x0,...,xn et la nullité des dérivées d’ordre (d−1)/2 à d/2 en xn. On résoud le système et on remplace pour avoir les valeurs numériques des coefficients du polynome.

4.7  Autres approximations polynomiales.

On peut citer le polynôme de Taylor en un point qui donne une bonne approximation près du point (voir aussi ci-dessus le polynôme de Hermite), les approximants de Padé où on approche par le quotient de 2 polynômes (ceci donne parfois de très bons résultats comme pour la fonction exponentielle près de 0 par exemple).

On observe que le dénominateur est le numérateur pris en −x. Ceci permet de calculer précisément ex par exemple sur [−1,1] en très peu d’opérations arithmétiques, calcul de X=xx, puis de a=(30X+3360)X+30240, b=x((X+420)X+15120) puis (a+b)/(ab), soit 5 multiplications, 6 additions, 1 division.

On utilise aussi souvent des approximations distinctes selon l’intervalle considéré, par exemple l’interpolation linéaire par morceaux, les fonctions splines ... ou la réduction d’arguments pour se ramener à un intervalle avec une bonne apprixomation (par exemple ex=(ex/2)2)

5  Intégration numérique

Les fractions rationnelles admettent une primitive que l’on calcule en décomposant la fraction avec Bézout comme expliqué précédemment. Mais elles font figure d’exceptions, la plupart des fonctions n’admettent pas de primitives qui s’expriment à l’aide des fonctions usuelles. Pour calculer une intégrale,on revient donc à la définition d’aire sous la courbe, aire que l’on approche, en utilisant par exemple un polynome de Lagrange.

Le principe est donc le suivant : on découpe l’intervalle d’intégration en subdivisions [a,b]=[a,a+h] + [a+h,a+2h]+...[a+(n−1)h,a+nh=b, où h=(ba)/n est le pas de la subdivision, et sur chaque subdivision, on approche l’aire sous la courbe.

5.1  Les rectangles et les trapèzes

Sur une subdivision [α,β], on approche la fonction par un segment. Pour les rectangles, il s’agit d’une horizontale : on peut prendre f(α), f(β) (rectangle à droite et gauche) ou f((α+β)/2) (point milieu), pour les trapèzes on utilise le segment reliant [α,f(α)] à [β,f(β)].

Exemple : calcul de la valeur approchée de ∫01 t3 dt (on en connait la valeur exacte 1/4=0.25) par ces méthodes en subdivisant [0,1] en 10 subdivisions (pas h=1/10), donc α=j/10 et β=(j+1)/10 pour j variant de 0 à 9. Pour les rectangles à gauche, on obtient sur une subdivision f(α)=(j/10)3 que l’on multiplie par la longueur de la subdivision soit h=1/10 :

1
10
 
9
j=0
 (
j
10
)3  = 
81
400
 = 0.2025

Pour les rectangles à droite, on obtient

1
10
 
10
j=1
 (
j
10
)3  = 
121
400
 = 0.3025

Pour le point milieu f((α+β)/2)=f((j/10+(j+1)/10)/2)=f(j/10+1/20)

1
10
 
9
j=0
 (
j
10
+
1
20
)3  = 199/800 = 0.24875

Enfin pour les trapèzes, l’aire du trapèze délimité par l’axe des x, les verticales y=α, y=β et les points sur ces verticales d’ordonnées respectives f(α) et f(β) vaut

h 
f(α)+f(β)
2

donc

1
10
 
9
j=0
 


(
j
10
)3 +(
j+1
10
)3


101
400
 = 0.2525 

Dans la somme des trapèzes, on voit que chaque terme apparait deux fois sauf le premier et le dernier.

Plus généralement, les formules sont donc les suivantes :

     
rectangle gauche =
 h 
n−1
j=0
 f(a+jh 
    (4)
rectangle droit =
 h 
n
j=1
 f(a+jh 
    (5)
point milieu =
 h 
n−1
j=0
 f(a+jh+
h
2
 
    (6)
trapezes  =
  h 


f(a)+f(b)
2
+
n−1
j=1
 f(a+jh


 
    (7)

h=(ba)/n est le pas de la subdivision, n le nombre de subdivisions.





On observe sur l’exemple que le point milieu et les trapèzes donnent une bien meilleure précision que les rectangles. Plus généralement, la précision de l’approximation n’est pas la même selon le choix de méthode. Ainsi pour les rectangles à gauche (le résultat est le même à droite), si f est continument dérivable, de dérivée majorée par une constante M1 sur [a,b], en faisant un développement de Taylor de f en α, on obtient

|
β


α
 f(tdt − 
β


α
 f(α) dt | = | 
β


α
 f′(θt)(t−α) dt | ≤ M1 
β


α
 (t−α) dt = M1
(β−α)2
2

Ainsi dans l’exemple, on a M1=3, l’erreur est donc majorée par 0.015 sur une subdivision, donc par 0.15 sur les 10 subdivisions.

Pour le point milieu, on fait le développement en (α+β)/2 à l’ordre 2, en supposant que f est deux fois continument dérivable :

 |
β


α
 f(t)  − 
β


α
 f(
α+β
2
)  |
=
β


α
 f′(
α+β
2
)(t
α+β
2
dt
  
+  
β


α
 
f′′(θt)
2
(t
α+β
2
)2 |
 
M2
2
 2 
β


α+β
2
(t
α+β
2
)2 dt  
 
M2
(β−α)3
24

Dans l’exemple, on a M2=6, donc l’erreur sur une subdivision est majorée par 0.25e−3, donc sur 10 subdivisions par 0.25e−2=0.0025.

Pour les trapèzes, la fonction g dont le graphe est le segment reliant [α,f(α)] à [β,f(β)] est f(α)+(t−α)/(β−α)f(β), c’est en fait un polynome de Lagrange, si f est deux fois continument dérivable, on peut donc majorer la différence entre f et g en utilisant (1), on intègre la valeur absolue ce qui donne

|
β


α
 f(tdt − 
β


α
 g(tdt | ≤ 
β


α
 |
f′′(ξx)
2
 (x−α)(x−β)| ≤ M2 
(β−α)3
12
 

M2 est un majorant de |f′′| sur [a,b].

Lorsqu’on calcule l’intégrale sur [a,b] par une de ces méthodes, on fait la somme sur n=(ba)/h subdivisions de longueur β−α=h, on obtient donc une majoration de l’erreur commise sur l’intégrale :

Lorsque h tend vers 0, l’erreur tend vers 0, mais pas à la même vitesse, plus rapidement pour les trapèzes et le point milieu que pour les rectangles. Plus on approche précisément la fonction sur une subdivision, plus la puissance de h va être grande, plus la convergence sera rapide lorsque h sera petit, avec toutefois une contrainte fixée par la valeur de Mk, borne sur la dérivée k-ième de f (plus k est grand, plus Mk est grand en général). Nous allons voir dans la suite comment se comporte cette puissance de h en fonction de la facon dont on approche f.

5.2  Ordre d’une méthode

On appelle méthode d’intégration l’écriture d’une approximation de l’intégrale sur une subdivision sous la forme

β


α
 f(tdt ≈ I(f)=
k
j=1
 wj f(xj

où les xj sont dans l’intervalle [α,β], par exemple équirépartis sur [α,β]. On utilise aussi la définition :

β


α
 f(tdt ≈ I(f)= (β−α)
k
j=1
 wj f(xj

On prend toujours ∑j wj=β−α (ou ∑j wj=1) pour que la méthode donne le résultat exact si la fonction est constante.

On dit qu’une méthode d’intégration est d’ordre n si il y a égalité ci-dessus pour tous les polynômes de degré inférieur ou égal à n et non égalité pour un polynôme de degré n+1. Par exemple, les rectangles à droite et gauche sont d’ordre 0, le point milieu et les trapèzes sont d’ordre 1. Plus généralement, si on approche f par son polynôme d’interpolation de Lagrange en n+1 points (donc par un polynôme de degré inférieur ou égal à n), on obtient une méthode d’intégration d’ordre au moins n.

Si une méthode est d’ordre n avec des wj≥ 0 et si f est n+1 fois continument dérivable, alors sur une subdivision, on a :

  |
β


α
 fI(f)| ≤ Mn+1 
(β−α)n+2
(n+1)!
(
1
n+2
+1)     (8)

En effet, on fait le développement de Taylor de f par exemple en α à l’ordre n

 f(t)=
Tn(f)+
(t−α)n+1
(n+1)!
 f[n+1]t),
 Tn(f)=
f(α)+(t−α)f′(α)+...+ 
(t−α)n
n!
 f[n](α)

Donc

|
β


α
 f− 
β


α
 Tn(f)| ≤ 
β


α
 
(t−α)n+1
(n+1)!
 |f[n+1]t)|  ≤ 


Mn+1 
(t−α)n+2
(n+2)!
 


β



α

De plus,

 |I(f) −I(Tn(f))| =|I


f[n+1]t)
(t−α)n+1
(n+1)!
 


|
k
j=1
 |wjMn+1 
(xj−α)n+1
(n+1)!
 
 
k
j=1
 |wjMn+1 
(β−α)n+1
(n+1)!

Donc comme la méthode est exacte pour Tn(f), on en déduit que

|
β


α
 fI(f)|
=
|
β


α
 f
β


α
 Tn(f)+I(Tn(f))− I(f)| 
 
|
β


α
 f
β


α
 Tn(f)|+|I(Tn(f))− I(f)|
 
Mn+1  
(β−α)n+2
(n+2)!
 +  
k
j=1
 |wjMn+1 
(β−α)n+1
(n+1)!
 

Si les wj≥ 0, alors ∑j=1k |wj|=∑j=1k wj=β−α et on obtient finalement (8)

On remarque qu’on peut améliorer la valeur de la constante en faisant tous les développement de Taylor en (α+β)/2 au lieu de α, Après sommation sur les n subdivisions, on obtient que :

Théorème 10   Pour une méthode d’ordre n à coefficients positifs et une fonction f n+1 fois continument dérivable 
|
b


a
 fI(f)| ≤ Mn+1 
hn+1
2n+1(n+1)!
  (ba)  (
1
(n+2)
+1)

On observe que cette majoration a la bonne puissance de h sur les exemples déja traités, mais pas forcément le meilleur coefficient possible, parce que nous avons traité le cas général d’une méthode d’ordre n, et utilisé une majoration pas toujours optimale du reste. Pour obtenir la meilleure valeur possible de la constante, il faut exprimer le reste de la formule de Taylor sous forme intégrale et utiliser la forme précise de la méthode :

f(t)=Tn(f)(t)+
t


α
 
(tx)n
n!
 f[n+1](xdx 

donc

β


α
 f− 
β


α
 Tn(f) = 
β


α
 
t


α
(tx)n
n!
  f[n+1](x)  dx  dt
β


α
 
β


α
(tx)+n
n!
  f[n+1](x)  dx  dt

où (tx)+=max(0,tx). En intervertissant les deux intégrales (Fubini), on obtient :

β


α
 f− 
β


α
 Tn(f)
=
β


α
 
β


α
(tx)+n
n!
   dt f[n+1](x)  dx 
 =
β


α
 
β


x
(tx)+n
n!
  dt  f[n+1](x)  dx
 =
β


α
 
(β−x)n+1
(n+1)!
 f[n+1](x)  dx

D’autre part :

I(f)−I(Tn(f))=
I




β


α
(xt)+n
n!
  f[n+1](t)  dt 




 =
k
j=1
 wj 
β


α
(xjt)+n
n!
 f[n+1](t)  dt 

Donc :

β


α
 fI(f) = 
β


α
  ( 
(β−x)n+1
(n+1)!
 −  
k
j=1
 wj 
(xjx)+n
n!
) f[n+1](x)  dx

On appelle noyau de Péano l’expression

N(x)=
(β−x)n+1
(n+1)!
 −  
k
j=1
 wj
(xjx)+n
n!

On a alors :

|
β


α
 fI(f)| ≤ 
β


α
 |N(x)| |f[n+1](x)|  dx ≤ Mn+1 
β


α
 |N(x)|  dx

5.3  Simpson

Il s’agit de la méthode obtenue en approchant la fonction sur la subdivision [α,β] par son polynome de Lagrange aux points α,(α+β)/2,β. On calcule l’intégrale par exemple avec un logiciel de calcul formel, avec Xcas :

qui donne la formule sur une subdivision

I(f) = 
h
6
 (f(α)+4f(
α+β
2
) + f(β)) 

et sur [a,b] :

  I(f) = 
h
6
 


f(a)+f(b)+ 4 
n−1
j=0
 f(a+jh+
h
2
) + 2  
n−1
j=1
 f(a+jh


    (9)



Si on intègre t3 sur [0,1] en 1 subdivision par cette méthode, on obtient

1
6
 (0+ 4 
1
23
 + 1)=
1
4
 

c’est-à-dire le résultat exact, ceci est aussi vérifié pour f polynome de degré inférieur ou égal à 2 puisque l’approximation de Lagrange de f est alors égale à f. On en déduit que la méthode de Simpson est d’ordre 3 (pas plus car la méthode de Simpson appliquée à l’intégrale de t4 sur [0,1] n’est pas exacte). On peut améliorer la constante générale de la section précédente pour la majoration de l’erreur en :

|
b


a
 f − I(f)| ≤ 
h4
2880
 (baM4 

En effet sur une subdivision élémentaire [α,β], le noyau de Péano vaut :

 N(x)=
(β−x)4
4!
 − 
1
6
 
(β−x)+3
3!
2
3
 
(
α+β
2
x)+3
3!
 − 
1
6
 
(α−x)+3
3!
 
 =
(β−x)4
4!
 − 
1
6
 
(β−x)3
3!
2
3
 
(
α+β
2
x)+3
3!
 

on observe que N(x)≤ 0 sur [α,β]




et son intégrale vaut −1/2880(β−α)4 :

La méthode de Simpson nécessite 2n+1 évaluations de f (le calcul de f est un point étant presque toujours l’opération la plus couteuse en temps d’une méthode de quadrature), au lieu de n pour les rectangles et le point milieu et n+1 pour les trapèzes. Mais on a une majoration en h4 au lieu de h2 donc le “rapport qualité-prix” de la méthode de Simpson est meilleur, on l’utilise donc plutot que les méthodes précédentes sauf si f n’a pas la régularité suffisante (ou si M4 est trop grand).

5.4  Newton-Cotes

On peut généraliser l’idée précédente, découper la subdivision [α,β] en n parts égales et utiliser le polynôme d’interpolation en ces n+1 points x0=α, x1, ..., xn=β. Ce sont les méthodes de Newton-Cotes, qui sont d’ordre n au moins. Comme le polynôme d’interpolation dépend linéairement des ordonnées, cette méthode est bien de la forme :

I(f)=(β−α)
n
j=0
 wj f(xj)

De plus les wj sont universels (ils ne dépendent pas de la subdivision), parce qu’on peut faire le changement de variables x=α+t(β−α) dans l’intégrale et le polynôme d’interpolation et donc se ramener à [0,1].

Exemple : on prend le polynôme d’interpolation en 5 points équidistribués sur une subdivision [a,b] (méthode de Boole). Pour calculer les wj, on se ramène à [0,1], puis on tape

et on lit les coefficients de f0 à f4 qui sont les w0 à w4: 7/90, 32/90, 12/90, 32/90, 7/90. Voir aussi la section 5.5 La méthode est d’ordre au moins 4 par construction, mais on vérifie qu’elle est en fait d’ordre 5 (exercice), la majoration de l’erreur d’une méthode d’ordre 5 est

|
b


a
 f −I(f)| ≤ 
M6
26 6!
(1+
1
7
h6 (ba

elle peut être améliorée pour cette méthode précise en

|
b


a
 f −I(f)| ≤ 
M6
1935360
 h6 (ba

En pratique, on ne les utilise pas très souvent, car d’une part pour n≥ 8, les wj ne sont pas tous positifs, et d’autre part, parce que la constante Mn devient trop grande. On préfère utiliser la méthode de Simpson en utilisant un pas plus petit.

Il existe aussi d’autres méthodes, par exemple les quadratures de Gauss (on choisit d’interpoler en utilisant des points non équirépartis tels que l’ordre de la méthode soit le plus grand possible, cf. infra) ou la méthode de Romberg qui est une méthode d’accélération de convergence basée sur la méthode des trapèzes (on prend la méthode des trapèzes en 1 subdivision de [a,b], puis 2, puis 22, ..., et on élimine les puissances de h du reste ∫fI(f) en utilisant un théorème d’Euler-Mac Laurin qui montre que le développement asymptotique de l’erreur en fonction de h ne contient que des puissances paires de h). De plus, on peut être amené à faire varier le pas h en fonction de la plus ou moins grande régularité de la fonction.

5.5  Calcul des poids wi

Si la méthode d’intégration consiste à interpoler f en n points x0,...,xn, alors la méthode est exacte pour tout polynôme de degré n. Si on prend Pj(x)=∏kj(xxk), on en déduit :

β


α
Pj(x)  dx = wj Pj(xj

Par exemple en interpolant en 0,1/2,1 sur [0,1], on obtient

w0 P0(0)=
1


0
 P0(x)  dx,    P0(x)=(x
1
2
)(x−1)


On peut aussi résoudre un système linéaire en prenant pour f les polynômes de la base canonique, la matrice du système est la transposée de la matrice de Vandermonde en les xj et le second membre a pour j-ième composante ∫αβxj dx.

5.6  En résumé

Intégration sur [a,b], h pas d’une subdivision, Mk majorant de la dérivée k-ième de la fonction sur [a,b]

 formuleLagrange degréordreerreur
rectangles(4), (5)00M1 h (ba)/2
point milieu(6)01M2 h2 (ba)/24
trapèzes(7)11M2 h2 (ba)/12
Simpson(9)23M4 h4 (ba)/2880

5.7  Accélération de Richardson-Romberg

Proposition 11   Soit g une fonction de classe C2k sur [a,b], Th(g) la valeur de la méthode des trapèzes sur [a,b] de pas h=(ba)/N (N entier). Alors Th(g) admet un développement en puissances paires de h à l’ordre 2k.

Pour montrer ce résultat, il faut établir la formule d’Euler-Mac Laurin. On commence par se placer sur une subdivision de l’intervalle [0,1], on intègre par parties ∫01 f(t) dt en faisant apparaitre la formule des trapèzes, on intègre donc 1 en t−1/2

1


0
 f(t)  dt =  [(t
1
2
)f(t)]01 − 
1


0
 (t
1
2
)f′(t)  dt =
f(0)+f(1)
2
 − 
1


0
 (t
1
2
)f′(t)  dt

où encore

T1(f)=
1


0
 f(t)  dt  +  
1


0
 (t
1
2
)f′(t)  dt 

Pour poursuivre, on pose B1(t)=t−1/2, qu’on va intégrer en 1/2((t−1/2)2+c), où on choisit c pour que l’intégrale soit nulle, donc c=−1/6. On pose B2=(t−1/2)2−1/6, on a :

T1(f)=
1


0
 f(t)  dt  + [
B2
2
 f′]01 − 
1


0
 
B2
2
f′′(t)  dt 

Plus généralement, on pose

Bp+1′=pBp,    
1


0
 Bp+1(t)  dt=0 

ce qui définit de manière unique les Bp. La nullité de l’intégrale montre que Bp+1(1)=Bp+1(0) ce qui simplifiera l’expression des termes tout intégrés. De plus, on montre par récurrence que les Bp ont une symétrie paire ou impaire selon la parité de p par rapport à t=1/2. Après p intégrations par parties, on obtient :

  T1(f)=
1


0
 f(t)  dt  + 
B2(0)
2
 (f′(1)−f′(0))  + ... 
  
B2k(0)
(2k)!
 (f[2k−1](1)−f[2k−1](0))  − 
1


0
 
B2k
(2k)!
 f[2k](t)  dt 

En faisant le même raisonnement sur [k,k+1] pour k=1,...,N−1 et en sommant, on obtient la formule d’Euler-Mac Laurin :

  T1[0,N](f)=
N


0
 f(t)  dt  + 
B2(0)
2
 (f′(N)−f′(0))  + ... 
  
B2k(0)
(2k)!
 (f[2k−1](N)−f[2k−1](0))  − 
N


0
 
B2k
(2k)!
 f[2k](t)  dt 

On pose alors x=a+ht (donc dt=dx/h) et f(t)=g(x) (donc f′(t)=df/dt=h dg/dx), on obtient

 
1
h
 Th[a,b](g)
=
1
h
 
N


0
 g(x)  dx  + 
B2(0)
2
 h(g′(b)−g′(a))  + ... 
  
B2k(0)
(2k)!
 h2k−1 (g[2k−1](N)−g[2k−1](0))  − 
N


0
 
B2k
(2k)!
 h2k g[2k](x)  
1
h
 dx 

donc

 Th[a,b](g)=
N


0
 g(x)  dx  + h2 
B2(0)
2
 (g′(b)−g′(a))  + ... 
  
h2k 
B2k(0)
(2k)!
  (g[2k−1](N)−g[2k−1](0))  − h2k 
N


0
 
B2k
(2k)!
  g[2k](x)   dx 

L’accélération consiste à éliminer les puissances de h2 en commençant par h2 avec des subdivisions deux fois plus fines à chaque itération. Ainsi Th1(f) = (4Th/2(f)−Th(f))/(4−1) n’a plus de termes en h2 et tend vers l’intégrale à approcher lorsque h tend vers 0. On peut d’ailleurs vérifier qu’il s’agit de la méthode de Simpson. On élimine ensuite le terme en h4 en posant Th2(f)=(42 Th/21(f)−Th1(f))/(42−1) et ainsi de suite. On construit un tableau triangulaire T dont chaque ligne l contient Th/2l(f), Th/2l1(f), ... (avec des indices qui commencent à 0). Pour calculer le terme d’indice 0 de la ligne courante on fait une méthode des trapèzes sur 2 fois plus de subdivisions que la précédente, puis pour le j-ième terme T[l,j] on effectue (4^j*T[l-1,j-1]-T[l,j-1])/(4^j-1) (on n’a donc besoin que de la ligne précédente pour calculer la ligne courante). On s’arrête par exemple lorsque la valeur absolue de la différence entre les derniers termes de deux lignes consécutives est inférieur à la précision souhaitée (erreur empirique).

5.8  Cas des fonctions périodiques.

Si f est une fonction périodique régulière (C), alors la méthode des trapèzes sur une période est d’ordre arbitrairement grand. En effet, pour une série de Fourier tronquée à l’ordre m, la formule des trapèzes avec N subdivisions donne le résultat exact de ∫0T f(t) dt dès que N>m. Il suffit ensuite d’utiliser que le reste de la série de Fourier (m>N) a des coefficients à décroissance rapide.

La méthode des trapèzes donne donc de bons résultats pour une fonction périodique, on peut d’ailleurs aussi l’appliquer pour calculer une valeur approchée des coefficients de Fourier de la fonction. La liste des valeurs approchées obtenue est alors la transformée de Fourier discrète des valeurs de la fonction f aux N points de la subdivision, elle se calcule donc rapidement avec la transformée de Fourier rapide. Par exemple, pour approcher f(x)=esin(x), on peut utiliser les commandes suivantes en Xcas :

f(x):=exp(sin(x));
N:=16; F:=seq(f(k/N*2.*pi),k,0,N-1); G:=fft(F);
k:=4;
g:=G[0]+sum(G[j]*exp(i*j*x),j,1,k)+sum(G[N-j]*exp(-i*j*x),j,1,k);
h:=normal(re(exp2trig(g)/N));



Ou directement 2*re(G[j]/N) est une valeur approchée du j-ième coefficient de Fourier aj de f, et -2*im(G[j]/N) de bj, par exemple :


On observe en effet une très bonne concordance (11 décimales).

Bien entendu, cela n’est pas très utile pour approcher esin(x) (il vaut mieux composer exponentielle et sinus), mais cela pourrait le devenir pour une fonction périodique plus compliquée ou pour une fonction périodique dont on ne connait qu’un échantillonage régulier (par exemple un fichier numérique audio).

5.9  Quadratures gaussiennes.

5.9.1  Description

On a vu que l’interpolation polynomiale était de meilleure qualité en prenant les points de Tchebyshev plutot que des points équidistants, il est donc naturel de calculer des approximations d’intégrale de cette manière ou encore d’optimiser le choix des abscisses pour avoir une méthode d’intégration d’ordre maximal.

Si on se fixe n abscisses x1 à xn, on peut obtenir l’ordre 2n−1. En effet, considérons le polynôme Pn=∏i=1n (xxi), qui est de degré n, si la méthode est d’ordre 2n−1 alors il sera orthogonal à tous les polynômes de degré inférieur strict à n pour le produit scalaire

<f|g>=
b


a
 f(xg(x)  dx

puisque <Pn|xj>=∫ab Pn xj sera combinaison linéaire des Pn xj en xk, k=1..n (car la formule d’intégration est exacte puisque le degré du polynôme Pnxj est au plus 2n−1). Donc Pn est à une constante multiplicative près le n-ième polynôme orthogonal pour l’integrale sur [a,b], si [a,b]=[−1,1] c’est legendre(n). Réciproquement, si les xk sont les racines de ce polynôme, alors la formule d’intégration est exacte, on effectue la division euclidienne du polynôme P de degré au plus 2n−1 à intégrer par Pn

PPn Q + R ,    deg(Q) ≤ n−1 

On a ∫ab Pn Q=0 par orthogonalité et la combinaison linéaire correspondante en les xk est nulle, et on a exactiture pour R, car de degré au plus n−1.

Exemple :



5.9.2  Calcul des poids

On peut calculer les poids en appliquant la section 5.5. On peut ainsi montrer que les poids sont positifs en appliquant la formule d’intégration au polynôme ∏jk(xxj) (la formule est exacte à cause du degré du polynome). On peut d’ailleurs montrer que le poids wi en une racine du polynôme de Legendre sur [−1,1] est donné par :

wi=
2
(1−xi2)Pn(xi)2
 = 
2(1−xi2)
n2Pn−1(xi)2
=
2
nPn−1(xi)Pn′(xi)






Preuve de la dernière formule :
On a I(Pn/(xxi))=wi Pn(xi) la valeur de (Pn/(xxi)) en xi est par définition de Pn(xi) la limite (Pn(x)−Pn(xi))/(xxi) lorsque x tend vers xi (rappelons que Pn(xi)=0). Par exactitude de la formule d’intégration

wi Pn(xi)=
1


−1
 
Pn
xxi
 =
1
Pn−1(xi)
1


−1
 Pn 
Pn−1(xi)
xxi
  = 
1
Pn−1(xi)
1


−1
 Pn 
Pn−1(x)
xxi
 

la dernière égalité résulte du fait que (Pn−1(x)−Pn−1(xi))/(xxi) est un polynôme de degré au plus n−2 donc orthogonal à Pn. Donc

wi Pn(xi)= = 
1
Pn−1(xi)
1


−1
 Pn−1 
Pn(x)
xxi
 

Or Pn(x)/(xxi) est un polynôme de degré n−1, si on fait son quotient par Pn−1 on obtient une constante qui vaut an−1 en utilisant la relation de récurrence à 2 crans Pn=(an−1xbn−1)Pn−1cn−1Pn−2. Le reste est de degré n−2 donc orthogonal à Pn−1, d’où

wi Pn(xi)= = 
1
Pn−1(xi)
1


−1
 Pn−12 an−1 

On conclut en utilisant la valeur de an−1=(2n−1)/n et de || Pn−1||2=2/(2(n−1)+1).

5.9.3  Erreur d’une quadrature gaussienne

On considère une quadrature gaussienne d’ordre 2n+1 obtenue par interpolation aux n+1 racines du n+1-ième polynome de Legendre sur [−1,1].

Rappelons que le polynôme d’interpolation de f aux points x0,x1,...,xn s’écrit

Pn(t)=f[x0]+f[x0,x1](tx0)+...+f[x0,...,xn](tx0)...(txn−1)

Si on ajoute x on obtient un polynôme d’interpolation de degré n+1

Pn+1(t)=f[x0]+f[x0,x1](tx0)+...+f[x0,...,xn,x](tx0)...(txn−1)(txn)

qui coincide avec f en t=x donc

f(x)=f[x0]+f[x0,x1](xx0)+...+f[x0,...,xn,x](xx0)...(xxn−1)(xxn)

Le même calcul fait avec un point xn+1 en plus donne

f(x)=f[x0]+f[x0,x1](xx0)+...+f[x0,...,xn,xn+1](xx0)...(xxn−1)(xxn) + f[x0,...,xn,xn+1,x](xx0)...(xxn−1)(xxn)(xxn+1)

On a donc montré la formule :

  f[x0,...,xn,x]=f[x0,...,xn,xn+1]+f[x0,...,xn,xn+1,x](xxn+1)     (10)

mais en fait on va ajouter plus qu’un point, on va en ajouter n+1

  f(x)−(f[x0]+f[x0,x1](xx0)+...+f[x0,...,xn](xx0)...(xxn−1) )
 =f[x0,...,xn,xn+1](xx0)...(xxn−1)(xxn) +...+f[x0,...,x2n+1](xx0)...(xx2n) + f[x0,...,x2n+1,x](xx0)...(xx2n+1)

puis on intègre sur [−1,1], on obtient à gauche l’erreur et à droite, seul le dernier terme est non nul, car le polynôme de Legendre proportionnel à (xx0)...(xxn−1)(xxn) , de degré n+1, est orthogonal à tous les polynôme de degré ≤ n

1


−1
 (f(x) −Pn(x))  dx = 
1


−1
f[x0,...,x2n+1,x](xx0)...(xx2n+1)  dx 

D’autre part, le résultat sur l’erreur d’interpolation 1 donne

f[x0,...,xn,x]=
f[n+1]x)
(n+1)!

en particulier |f[x0,...,xn,x]| ≤ Mn+1/(n+1)! D’où la majoration

|  
1


−1
 (f(x) −Pn(x))  dx | ≤ 
M2n+2
(2n+2)!
1


−1
|(xx0)...(xx2n+1)|  dx 

Il suffit ensuite de faire tendre les xn+1,...,x2n+1 vers x0,...,xn pour enlever la valeur absolue et obtenir

|  
1


−1
 (f(x) −Pn(x))  dx | ≤ 
M2n+2
(2n+2)!
1


−1
 Ln+12(x)  dx 

Ln+1 est proportionnel au n+1-ième polynôme de Legendre de coefficient dominant 1.

Par exemple, pour n=1, on a 2 points d’interpolation en les racines de L2=x2−1/3 et l’erreur d’interpolation est majorée par M4/24∫−11(x2−1/3) dx=M4/24 × 8/45=M4/135.

5.10  Méthode adaptative.

On calcule une valeur approchée de l’intégrale sur [a,b] par deux quadratures gaussiennes emboitées, on estime l’erreur, si elle est supérieure à la tolérance on divise en 2. On recommence en subdivisant en 2 l’intervalle où l’erreur est maximale. On s’arrête lorsque l’erreur estimée est inférieure à la tolérance.

L’estimation de l’erreur se fait par exemple avec deux quadratures gaussiennes emboitées (c’est-à-dire que les points d’interpolation de la moins fine sont contenues dans les points d’interpolation de la plus fine, pour éviter de devoir calculer la fonction en de nouveaux points, on considére alors l’erreur sur la quadrature la moins fine comme la valeur absolue de la différence des deux valeurs). Ou avec trois quadratures emboitées, Hairer propose de prendre comme quadrature la plus fine d’ordre 30 (15 points), intermédiaire d’ordre 14 (avec les mêmes points sauf le point central), moins fine d’ordre 6 (avec les points 1, 3, 5, 9, 11, 13), et d’estimer l’erreur par

є1=| I30I14|,  є2=| I30I6|;  є = є1 


є1
є2
 


2



 
 

On observe en effet que є est en h30, comme l’ordre de la méthode.

5.11  Accélération de Richardson-Romberg

Proposition 12   Soit g une fonction de classe C2k sur [a,b], Th(g) la valeur de la méthode des trapèzes sur [a,b] de pas h=(ba)/N (N entier). Alors Th(g) admet un développement en puissances paires de h à l’ordre 2k.

Pour montrer ce résultat, il faut établir la formule d’Euler-Mac Laurin. On commence par se placer sur une subdivision de l’intervalle [0,1], on intègre par parties ∫01 f(t) dt en faisant apparaitre la formule des trapèzes, on intègre donc 1 en t−1/2

1


0
 f(t)  dt =  [(t
1
2
)f(t)]01 − 
1


0
 (t
1
2
)f′(t)  dt =
f(0)+f(1)
2
 − 
1


0
 (t
1
2
)f′(t)  dt

où encore

T1(f)=
1


0
 f(t)  dt  +  
1


0
 (t
1
2
)f′(t)  dt 

Pour poursuivre, on pose B1(t)=t−1/2, qu’on va intégrer en 1/2((t−1/2)2+c), où on choisit c pour que l’intégrale soit nulle, donc c=−1/6. On pose B2=(t−1/2)2−1/6, on a :

T1(f)=
1


0
 f(t)  dt  + [
B2
2
 f′]01 − 
1


0
 
B2
2
f′′(t)  dt 

Plus généralement, on pose

Bp+1′=pBp,    
1


0
 Bp+1(t)  dt=0 

ce qui définit de manière unique les Bp. La nullité de l’intégrale montre que Bp+1(1)=Bp+1(0) ce qui simplifiera l’expression des termes tout intégrés. De plus, on montre par récurrence que les Bp ont une symétrie paire ou impaire selon la parité de p par rapport à t=1/2. Après p intégrations par parties, on obtient :

  T1(f)=
1


0
 f(t)  dt  + 
B2(0)
2
 (f′(1)−f′(0))  + ... 
  
B2k(0)
(2k)!
 (f[2k−1](1)−f[2k−1](0))  − 
1


0
 
B2k
(2k)!
 f[2k](t)  dt 

En faisant le même raisonnement sur [k,k+1] pour k=1,...,N−1 et en sommant, on obtient la formule d’Euler-Mac Laurin :

  T1[0,N](f)=
N


0
 f(t)  dt  + 
B2(0)
2
 (f′(N)−f′(0))  + ... 
  
B2k(0)
(2k)!
 (f[2k−1](N)−f[2k−1](0))  − 
N


0
 
B2k
(2k)!
 f[2k](t)  dt 

On pose alors x=a+ht (donc dt=dx/h) et f(t)=g(x) (donc f′(t)=df/dt=h dg/dx), on obtient

 
1
h
 Th[a,b](g)
=
1
h
 
N


0
 g(x)  dx  + 
B2(0)
2
 h(g′(b)−g′(a))  + ... 
  
B2k(0)
(2k)!
 h2k−1 (g[2k−1](N)−g[2k−1](0))  − 
N


0
 
B2k
(2k)!
 h2k g[2k](x)  
1
h
 dx 

donc

 Th[a,b](g)=
N


0
 g(x)  dx  + h2 
B2(0)
2
 (g′(b)−g′(a))  + ... 
  
h2k 
B2k(0)
(2k)!
  (g[2k−1](N)−g[2k−1](0))  − h2k 
N


0
 
B2k
(2k)!
  g[2k](x)   dx 

L’accélération consiste à éliminer les puissances de h2 en commençant par h2 avec des subdivisions deux fois plus fines à chaque itération. Ainsi Th1(f) = (4Th/2(f)−Th(f))/(4−1) n’a plus de termes en h2 et tend vers l’intégrale à approcher lorsque h tend vers 0. On peut d’ailleurs vérifier qu’il s’agit de la méthode de Simpson. On élimine ensuite le terme en h4 en posant Th2(f)=(42 Th/21(f)−Th1(f))/(42−1) et ainsi de suite. On construit un tableau triangulaire T dont chaque ligne l contient Th/2l(f), Th/2l1(f), ... (avec des indices qui commencent à 0). Pour calculer le terme d’indice 0 de la ligne courante on fait une méthode des trapèzes sur 2 fois plus de subdivisions que la précédente, puis pour le j-ième terme T[l,j] on effectue (4^j*T[l-1,j-1]-T[l,j-1])/(4^j-1) (on n’a donc besoin que de la ligne précédente pour calculer la ligne courante). On s’arrête par exemple lorsque la valeur absolue de la différence entre les derniers termes de deux lignes consécutives est inférieur à la précision souhaitée (erreur empirique).

5.12  Méthodes probabilistes.

Pour déterminer ∫ab f(t) dt , on l’interprète comme une espérance, plus précisément comme (ba)E(f(X)) où X est une variable aléatoire qui suit la loi uniforme sur [a,b], et on approche cette valeur par

ba
n
 
n
k=1
 f(xk

xk est obtenu par un générateur pseudo-aléatoire (selon la loi uniforme). Par exemple



La convergence en fonction de n est assez lente, on peut l’observer en faisant plusieurs estimations :



En effet, les tirages sont équidistribués selon la même loi, la loi des grands nombres s’applique donc : on fait ba fois une moyenne de n tirages, si n est grand, on converge vers une loi normale dont l’écart-type est en (ba)σ/√n. La valeur de la constante σ peut se calculer à partir de f

σ2 = 
1
ba
b


a
 f2(t)  dt −




1
ba
b


a
 f(t)  dt  




2





 
 

par exemple ici


mais on ne fait pas ce calcul en pratique (puisqu’il faudrait calculer une intégrale), on estime l’écart-type σ/√n de la loi normale par l’écart-type de l’échantillon des estimations stddevp(I).


On peut donc obtenir rapidement une estimation de σ en prenant l’écart-type d’une séquence de valeurs de f

Cette méthode converge donc beaucoup moins vite que les quadratures, en dimension 1. Mais elle se généralise très facilement en dimension plus grande en conservant la même vitesse de convergence alors que le travail nécessaire pour une méthode de quadrature croit comme une puissance de la dimension, et ne nécessite pas de paramétrer des domaines d’intégration compliqués (il suffit par exemple d’utiliser la méthode du rejet pour avoir un générateur uniforme dans un domaine inclus dans un cube).

6  Suites itératives et applications

Résumé:
Théorème du point fixe, méthode de Newton,convexité. Exemple: calcul de valeur approchée de racines carrées, Résolution d’équations.

6.1  Le point fixe dans ℝ

Soit f une fonction continue sur un intervalle I=[a,b] de ℝ, et à valeurs dans I (attention à bien choisir I pour que l’image de I par f reste dans I). On s’intéresse à la suite

  un+1=f(un),    u0 ∈ I      (11)

Supposons que un converge vers une limite lI lorsque n → +∞, alors la limite doit vérifier

f(l)=l 

puisque f est continue. On dit que l est un point fixe de f. Ceci amène à l’idée d’utiliser ces suites pour résoudre numériquement l’équation f(x)=x. Nous allons donner un théorème permettant d’assurer que la suite (11) converge, et que la limite est l’unique solution de f(l)=l sur I.

Définition 3   On dit que f est contractante de rapport k<1 sur I si
∀ x,y ∈ I,    |f(y)−f(x)| ≤ k |yx

En pratique, les fonctions f que l’on considèrera seront continument dérivables, donc d’après le théorème des accroissements finis

f(y)−f(x)=f′(θ) (yx),    θ ∈ [x,y

ainsi pour vérifier que f est contractante, on étudie la valeur absolue de f′ sur I, il suffit de montrer que cette valeur absolue est strictement inférieure à un réel k<1 pour conclure (il faut donc chercher le maximum de |f′| sur I. Attention, il s’agit du maximum de |f′| et pas du maximum de f′, ce qui revient à chercher le maximum de f′ et de −f′).

On a alors le

Théorème 13   (du point fixe)
si
f est contractante de I=[a,b] dans I de rapport k alors la suite (11) converge vers l’unique solution de f(l)=l dans I. On a de plus les encadrements :
   |unl| ≤ kn |ba|,    |un −l | ≤ 
|un+1un|
1−k
      (12)

Démonstration : Tout d’abord si f est contractante, on montre à partir de la définition de la continuité que f est continue. Soit g(x)=f(x)−x, alors g est continue, positive en a et négative en b, il existe donc l∈[a,b] tel que g(l)=0 (théorème des valeurs intermédiaires). Soit un une suite définie par (11). On a alors pour tout n

|un+1l|=|f(un)−f(l)| ≤ k |unl

Donc par une récurrence évidente :

|unl| ≤ kn |u0l

ce qui entraine d’ailleurs que |unl| ≤ kn |ab|. Comme k ∈ [0,1[ , la suite géométrique kn converge vers 0 lorsque n tend vers l’infini, donc un tend vers l. Notons que l est unique car si l′ est une autre solution alors |ll′|=|f(l)−f(l′)| ≤ k|ll′| donc (1−k)|ll′| ≤ 0, or 1−k>0 et |ll′| ≥ 0 donc |ll′| doit être nul. Passons à la preuve de la majoration (12) qui est importante en pratique car elle donne un test d’arrêt de calcul des termes de la suite récurrente, on écrit pour m>0 :

unlun − un+1 + un+1 − un+2 + ... + un+m−1− un+muml 

puis on majore avec l’inégalité triangulaire

|unl| ≤ 
m−1
j=0
 |un+jun+j+1| + |uml

puis on applique le fait que f est contractante de rapport k

|unl| ≤ 
m−1
j=0
 kj |unun+1| + |uml

soit

|unl| ≤ 
1−km
1−k
 |unun+1| + |uml

On fait alors tendre m vers l’infini d’où le résultat.

Exemples : Cherchons une valeur approchée de √2 par cette méthode. Il faut d’abord trouver une fonction f dont √2 est un point fixe, par exemple

f(x)=
x+2
x+1

On vérifie que f(√2)=√2), puis que f′=−1/(x+1)2 donc f décroit. On va voir si les hypothèses du théorème du point fixe s’appliquent sur par exemple [1,2]. Comme f est décroissante f([1,2])=[f(2),f(1)]=[4/3,3/2] qui est bien inclus dans [1,2] . De plus f′ est comprise entre −1/(1+1)2=−1/4 et −1/(2+1)2=−1/9 donc |f′|<1/4, f est contractante de rapport 1/4. On peut donc itérer la suite à partir par exemple de u0=1 et on va converger vers √2 (en s’en rapprochant à chaque cran d’un rapport inférieur à 1/4).

Considérons l’équation en x

x− e sin(x) =t,    e ∈ [0,1[ 

c’est l’équation du temps utilisée en astronomie pour trouver la position d’une planète sur son orbite elliptique (e étant l’excentricité de l’ellipse). Il n’y a pas de formule exacte permettant de calculer x en fonction de t. Si on a une valeur numérique pour t, on peut trouver une valeur numérique approchée de x par la méthode du point fixe, en réécrivant l’équation sous la forme

f(x)=t+esin(x) = x 

On observe que f envoie ℝ dans [te,t+e] donc on peut prendre I=[te,t+e], de plus |f′|≤ e <1, f est contractante de rapport e ∈ [0,1[, le théorème s’applique, il suffit de prendre une valeur initiale dans [te,t+e] et d’itérer la suite jusqu’à obtenir la précision désirée. Par exemple si on veut une valeur approchée de x à 10−6 près, il suffira que la différence entre deux termes successifs de la suite un vérifie

|un+1un| ≤  10−6 (1−e

on aura alors bien :

|unx| ≤ 
|un+1un|
1−e
 ≤ 10−6 

Cette méthode n’est pas toujours optimale, car la vitesse de convergence vers la limite l est dite “linéaire”, c’est-à-dire que le temps de calcul pour avoir n décimales est proportionnel à n (ou encore il faut effectuer un nombre d’itérations proportionnel à n, chaque itération faisant gagner en précision de l’ordre du rapport k de contractance). En effet, supposons que f′ est continue en l et que 0<L=|f′(l)|<1 . Il existe alors un intervalle I=[l−η,l+η] tel que

x ∈ I ⇒ 
L
2
 ≤ |f′(x)| ≤ 
1+L
2
 

Le théorème des accroissements finis donne alors

|un+1 − l | = |f(un)−f(l)| = |f′(θ)| |unl|,    θ ∈ [un,l

Si u0I, alors θ ∈ I donc |u1l| ≤ |u0l| et u1I, par récurrence on a pour tout n, unI

L
2
 |unl| ≤ |un+1 − l| ≤ 
1+L
2
 |unl

on a donc par récurrence




L
2



n



 
|u0l| ≤ |unl| ≤  


1+L
2
 


n



 
|u0l|

Donc pour avoir |unl| ≤ є il suffit que




1+L
2
 


n



 
|u0l| ≤ є ⇒ n ≥ 
ln(
є
|u0l|
)
ln( 
1+L
2
 

et il faut que




L
2



n



 
 |u0l| ≤ є ⇒ n ≥ 
ln(
є
|u0l|
)
ln( 
L
2

Si f est suffisamment régulière, il existe une méthode plus rapide lorsqu’on est proche de la racine ou lorsque la fonction a des propriétés de convexité, c’est la méthode de Newton. Et même si Newton n’est pas applicable, une simple dichotomie peut être plus efficace si la constante de contractance est supérieure à 1/2 (y compris prés de la solution de f(x)=x). Toutefois la méthode du point fixe reste intéressante si la constante de contractance est suffisamment petite (par exemple k=0.1 garantit 15 décimales en 15 itérations) et présente l’avantage de se généraliser en dimension plus grande, cf. la section suivante.

6.2  Le point fixe dans ℝn

Le théorème précédent se généralise.

Théorème 14   Soit I un ensemble fermé de n (ou d’un espace métrique complet) tel que f envoie I dans I et tel que f soit contractante sur I
∃ k<1, ∀ x,y ∈ I,    |f(x)−f(y)| ≤ k |xy|
Alors pour tout u0I, la suite (un) définie par un+1=f(un) converge vers l’unique solution dans I de f(l)=l.

La démonstration de la convergence est un peu différente de celle donnée en dimension 1, on montre que (un) est une suite de Cauchy, car pour n>m

|unum| ≤ 
n−1
j=m
|uj+1uj |  ≤ km 
1
1−k
 |u1u0

donc (un) est convergente puisque nous sommes dans un fermé d’un espace complet. (Cela permet d’ailleurs de généraliser l’énoncé donné en dimension 1 au cas où a ou b est infini)

La vitesse de convergence est linéaire, la démonstration est identique à celle de la dimension 1.

Remarque :

6.3  La méthode de Newton dans ℝ.

La méthode de Newton est une méthode de résolution de l’équation f(x)=0, attention à la différence avec le théorème du point fixe qui permet de résoudre numériquement f(x)=x. Si x0 est proche de la racine r on peut faire un développement de Taylor à l’ordre 1 de la fonction f en x0 :

f(x)=f(x0)+(xx0)f′(x0)+O((xx0)2

Pour trouver une valeur approchée de r, on ne garde que la partie linéaire du développement, on résout :

f(r)=0 ≈ f(x0) + (rx0f′(x0

donc (si f′(x0)≠ 0) :

r ≈ x0 −
f(x0)
f′(x0)

Graphiquement, cela revient à tracer la tangente à la courbe représentative de f et à chercher où elle coupe l’axe des x. On considère donc la suite récurrente définie par une valeur u0 proche de la racine et par la relation :

un+1 = un −
f(un)
f′(un)

Il y a deux théorèmes importants, l’un d’eux prouve que si u0 est “assez proche” de r alors la suite un converge vers r, malheureusement il est difficile de savoir en pratique si on est “assez proche” de u0 pour que ce théorème s’applique. Le second théorème donne un critère pratique facile à vérifier qui assure la convergence, il utilise les propriétés de convexité de la fonction.

Théorème 15   Soit f une fonction de classe C2 (2 fois continument dérivable) sur un intervalle fermé I. Soit r une racine simple de f située à l’intérieur de I (telle que f(r)=0 et f′(r)≠ 0). Alors il existe ε>0 tel que la suite définie par
un+1 = un −
f(un)
f′(un)
,    |u0r| ≤ ε 
converge vers r.

Si on a |f′′| ≤ M et |1/f′| ≤ m sur un intervalle [r−η,r+η] contenu dans I, alors on peut prendre tout réel ε>0 tel que ε < 2/(mM) et ε ≤ η.

Démonstration : on a

un+1r = un − r − 
f(un)
f′(un)
 = 
(unr)f′(un)−f(un)
f′(un)
 

En appliquant un développement de Taylor de f en un à l’ordre 2, on obtient pour un réel θ situé entre r et un :

0 = f(r)=f(un)+(runf′(un) + (run)2 
f′′(θ)
2
 

donc :

(unr)f′(un)−f(un)= (unr)2 
f′′(θ)
2
 

d’où :

|un+1r| ≤ |unr|2 
1
|f′(un)|
 
|f′′(θ)|
2
 

On commence par choisir un intervalle [r−ε,r+ε] contenant strictement r et tel que |f′′|<M et |1/f′|<m sur [r−ε,r+ε] (c’est toujours possible car f′′ et 1/f′ sont continues au voisinage de r puisque f′(r)≠ 0). Si un est dans cet intervalle, alors θ aussi donc

|un+1r| ≤ |unr|2 
Mm
2
 ≤  
|unr|Mm
2
|unr|,       (13)

On a |unr| ≤ ε, on diminue si nécessaire ε pour avoir ε < 2/(Mm), on a alors :

|un+1r| ≤ k |unr|,    k=
ε Mm
2
<1  

donc d’une part un+1 est encore dans l’intervalle [r−ε,r+ε] ce qui permettra de refaire le même raisonnement au rang suivant, et d’autre part on a une convergence au moins géométrique vers r. En fait la convergence est bien meilleure lorsqu’on est proche de r grace au carré dans |unr|2, plus précisément, on montre par récurrence que

|unr| ≤ |u0 − r|2n 


Mm
2
 


2n−1



 

il faut donc un nombre d’itérations proportionnel à ln(n) pour atteindre une précision donnée.

Remarque : ce théorème se généralise sur ℂ et même sur ℝn (cf. la section suivante).

Exemple : pour calculer √2, on écrit l’équation x2−2=0 qui a √2 comme racine simple sur I=[1/2,2], on obtient la suite récurrente

un+1 = un − 
un2−2
2un
 

Si on prend η=1/2, on a f′=2x et f′′=2 donc on peut prendre M=2 et m=1 car |1/f′|≤ 1 sur [√2−1/2,√2+1/2]. On a 2/(mM)=1, on peut donc prendre ε=1/2, la suite convergera pour tout u0 ∈ [√2−1/2,√2+1/2].

Plus généralement, on peut calculer une racine k-ième d’un réel a en résolvant f(x)=xka par la méthode de Newton.

L’inconvénient de ce théorème est qu’il est difficile de savoir si la valeur de départ qu’on a choisie se trouve suffisamment près d’une racine pour que la suite converge. Pour illustrer le phénomène, on peut par exemple colorer les points du plan complexe en n+1 couleurs selon que la suite définie par la méthode de Newton converge vers l’une des n racines d’un polynôme de degré n fixé au bout de par exemple 50 itérations (la n+1-ième couleur servant aux origines de suite qui ne semblent pas converger).

Passons maintenant à un critère très utile en pratique :

Définition 4   (convexité)
Une fonction
f continument dérivable sur un intervalle I de est dite convexe si son graphe est au-dessus de la tangente en tout point de I.

Il existe un critère simple permettant de savoir si une fonction de classe C2 est convexe :

Théorème 16   Si f est C2 et f′′ ≥ 0 sur I alors f est convexe.

Démonstration :
L’équation de la tangente au graphe en x0 est

y=f(x0)+f′(x0)(xx0

Soit

g(x)=f(x)−(f(x0)+f′(x0)(xx0)) 

on a :

g(x0)=0,    g′(x)=f′(x)−f′(x0),    g′(x0)=0,     g′′=f′′ ≥ 0 

donc g′ est croissante, comme g′(x0)=0, g′ est négative pour x<x0 et positive pour x>x0, donc g est décroissante pour x<x0 et croissante pour x>x0. On conclut alors que g ≥ 0 puisque g(x0)=0. Donc f est bien au-dessus de sa tangente.

On arrive au deuxième théorème sur la méthode de Newton

Théorème 17   Si f(r)=0, f′(r)>0 et si f′′ ≥ 0 sur [r,b] alors pour tout u0 ∈ [r,b] la suite de la méthode de Newton
un+1 = un −
f(un)
f′(un)
,  
est définie, décroissante, minorée par r et converge vers r. De plus
0 ≤ un −r ≤ 
f(un)
f′(r)
 

Démonstration :
On a f′′ ≥ 0 donc si f′(r)>0 alors f′>0 sur [r,b], f est donc strictement croissante sur [r,b] on en déduit que f>0 sur ]r,b] donc un+1un. Comme la courbe représentative de f est au-dessus de la tangente, on a un+1r (car un+1 est l’abscisse du point d’intersection de la tangente avec l’axe des x). La suite un est donc décroissante minorée par r, donc convergente vers une limite lr. À la limite, on a

l=l
f(l)
f′(l)
 ⇒ f(l)=0 

donc l=r car f>0 sur ]r,b].

Comme (un) est décroissante, on a bien 0 ≤ unr, pour montrer l’autre inégalité, on applique le théorème des accroissements finis, il existe θ ∈ [r,un] tel que

f(un)−f(r)=(unr)f′(θ) 

comme f(r)=0, on a

unr = 
f(un)
f′(θ)
 

et la deuxième inégalité du théorème en découle parce que f′ est croissante.

Variantes :
Il existe des variantes, par exemple si f′(r)<0 et f′′ ≥ 0 sur [a,r]. Si f′′ ≤ 0, on considère g=−f.

Application :
On peut calculer la valeur approchée de la racine k-ième d’un réel a>0 en appliquant ce deuxième théorème. En effet si a>0, alors xka est 2 fois continument dérivable et de dérivée première kxk−1 et seconde k(k−1)xk−2 strictement positives sur ℝ+∗ (car k ≥ 2). Il suffit donc de prendre une valeur de départ u0 plus grande que la racine k-ième, par exemple 1+a/k (en effet (1+a/k)k ≥ 1+k a/k=1+a). En appliquant l’inégalité du théorème, on a :

0 ≤ un − 

a
1
k
 
 ≤  
unk − a
k

a
1
k
 
k−1 
≤ 
unka
ka
  

a
1
k
 
≤ 
unka
ka
 (1+
a
k
)

Pour avoir une valeur approchée de (a)1/k à ε près, on peut donc choisir comme test d’arrêt

unk −a ≤ 
ka
1+
a
k
 ε 

Par exemple pour √2, le test d’arrêt serait un2−2 ≤ 2 ε.

6.4  La méthode de Newton dans ℝn.

Le premier énoncé du cas de la dimension 1 se généralise en :

Théorème 18   Soit f une fonction de classe C2 (2 fois continument dérivable) sur un fermé I de n. Soit r une racine simple de f située à l’intérieur de I (telle que f(r)=0 et f′(r)=(∂j fi)(r) inversible). Alors il existe ε>0 tel que la suite définie par
un+1 = un −(f′(un) )−1f(un),    |u0r| ≤ ε 
converge vers r.

Si on a |f′′| ≤ M et |(f′)−1| ≤ m sur une boule centré en r de rayon η>0 contenue dans I, alors on peut prendre tout réel ε>0 tel que ε < 2/(mM) et ε ≤ η.

La démonstration est calquée sur la dimension 1, mais il faut prendre le reste intégral dans la formule de Taylor

un+1r = un − r − f′(un)−1 f(un) =  f′(un)−1 (f′(un)(unr)−f(un) ) 

puis on applique Taylor le long du segment [r,un] :

0 = f(r)=f(un)+ f′(un) (run)+  
1


0
 (1−θ)(run) f′′(r+θ (unr)) (run)   dθ 

donc :

un+1r =−f′(un)−1  (run




1


0
 (1−θ)f′′(r+θ (unr))   dθ 




(run)  

et on en déduit (13) et on conclut de même en remplaçant intervalle centré en r de rayon ε par boule de rayon ε.

Remarque : la convergence “numérique” (au sens du calcul en flottant) de la suite un ne suffit pas à montrer l’existence d’une racine proche de un. Une méthode de preuve alternative au calcul des constantes m et M consiste à trouver un rectangle ou une boule autour de un préservée par l’application xxf′(x)−1 f(x).

6.5  Calcul approché des racines complexes simples

La section précédente nous a montré qu’on pouvait se ramener à la recherche de racines simples, ce qui donne envie d’essayer la méthode de Newton. On a malheureusement rarement la possibilité de pouvoir démontrer qu’à partir d’une valeur initiale donnée, la méthode de Newton converge, parce que les racines peuvent être complexes, et même si elles sont réelles, on n’a pas forcément de résultat sur la convexité du polynôme (cf. cependant une application des suites de Sturm qui permet de connaitre le signe de P′′ sur un intervalle sans le factoriser).

On effectue donc souvent des itérations de Newton, en partant de 0.0, en espérant s’approcher suffisamment d’une racine pour que le théorème de convergence théorique s’applique. On se fixe un nombre maximal d’itérations, si on le dépasse on prend alors une valeur initiale aléatoire complexe et on recommence.

Une fois une racine déterminée, on l’élimine en calculant le quotient euclidien Q de P par Xr (par l’algorithme de Horner), puis on calcule les racines du quotient Q (qui sont des racines de P).

Un problème pratique apparait alors, c’est que r n’est pas exact donc le quotient Q non plus, au fur et à mesure du calcul des racines de P, on perd de plus en plus de précision. Il existe une amélioration simple, si r′ est une racine approchée de Q, alors elle est racine approchée de P et on a toutes les chances qu’elle soit suffisamment proche d’une racine de P pour que le théorème s’applique, on effectue alors 1 ou 2 itérations de Newton avec r′ mais pour P (et non Q) afin d’améliorer sa précision comme racine de P.

Une méthode de calcul plus stable utilise la recherche des valeurs propres de la matrice companion en double précision, puis affine par la méthode de Newton pour obtenir des valeurs approchées multi-précision, c’est ce que fait proot, par exemple proot(x^3+x+1,50).

Enfin, on peut appliquer directement la méthode de Newton pour trouver dans ℂn toutes les racines simultanément, c’est la méthode de Durand-Kerner, Weierstrass. On pose gx(z)=∏i=1n (xzi), il s’agit de résoudre en z gx(z)=P(x). On a a l’ordre 1 en z

gx(z+w)=gx(z) − 
n
i=1
 wi 
 
j ≠ i
 (xzj) +O(w2) =P(x)

pour trouver wi, on pose x=zi, on obtient

 
j ≠ i
(zizjwi = P(zi)

donc

wi=−
P(zi)
 
j ≠ i
(zizj)

On peut aussi calculer le produit du dénominateur en effectuant gx′(zi) (la dérivée porte sur x). On retrouve la méthode de Newton à une variable où la dérivée du polynôme au dénominateur est remplacée par la valeur approchée du polynôme. D’où le programme

dw(P,N,eps):={ // Weierstrass, Durand-Kerner polynomial rooter
  local l,v,w,n,j,k,q,q1;
  P:=P/lcoeff(P);
  n:=degree(P);
  assume(l,symbol);
  v:=seq(exp(i*l/n*2.0*pi),l,0,n-1); w:=v;
  for k from 1 to N do
    q:=pcoeff(v);
    q1:=q';
    for j from 0 to n-1 do
      w[j]:=v[j]-horner(P,v[j])/horner(q1,v[j]);
    od;
    if (l2norm(w-v)<eps*l2norm(v))
      return w;
    v:=w;
  od;
  retourne "max iter reached";
}:;

Par exemple dw(x^3+x+1,100,1e-10) renvoie des valeurs approchées des racines de x3+x+1.

Si on s’intéresse seulement à la racine de module maximal d’un polynôme, on peut en trouver une estimation assez simplement en appliquant la méthode de la puissance à la matrice companion du polynôme. On peut améliorer la précision d’une racine par des itérations inverses ou par la méthode de Newton en une variable.

7  Réduction approchée des endomorphismes

On pourrait trouver des valeurs propres approchées d’une matrice en calculant le polynome caractéristique ou minimal puis en le factorisant numériquement. Mais cette méthode n’est pas idéale relativement aux erreurs d’arrondis (calcul du polynome caractéristiaue, de ses racines, et nouvelle approximation en calculant le noyau de A−λ I), lorsqu’on veut calculer quelques valeurs propres on préfère utiliser des méthodes itératives directement sur A ce qui évite la propagation des erreurs d’arrondi.

7.1  Méthode de la puissance

Elle permet de déterminer la plus grande valeur propre en valeur absolue d’une matrice diagonalisable lorsque celle-ci est unique. Supposons en effet que les valeurs propres de A soient x1,...,xn avec |x1| ≤ |x2| ≤ ... ≤ |xn−1| < |xn| et soient e1,...,en une base de vecteurs propres correspondants. On choisit un vecteur aléatoire v et on calcule la suite vk=Avk−1=Ak v . Si v a pour coordonnées V1,...,Vn) dans la base propre, alors

vk = 
n
j=1
 Vj xjk ej  = xnk wk,    wk=
n
j=1
 Vj 


xj
xn



k



 
 ej

L’hypothèse que xn est l’unique valeur propre de module maximal entraine alors que limk → +∞ wk = Vn en puisque la suite géométrique de raison xj/xn converge vers 0. Autrement dit, si Vn≠ 0 (ce qui a une probabilité 1 d’être vrai pour un vecteur aléatoire), vk est équivalent à Vn xnk en. Lorsque n est grand, vk est presque colinéaire au vecteur propre en (que l’on peut estimer par vk divisé par sa norme), ce que l’on détecte en testant si vk+1 et vk sont presques colinéaires. De plus le facteur de colinéarité entre vk+1 et vk est presque xn, la valeur propre de module maximal.

Exercice : tester la convergence de vk=Akv vers l’espace propre associé à λ=3 pour la matrice [[1,-1],[2,4]] et le vecteur v=(1,0). Attention à ne pas calculer Ak pour déterminer vk, utiliser la relation de récurrence!

Si on n’observe pas de convergence ou si elle est trop lente, alors |xn−1| est proche de |xn| ou égal, il est judicieux de faire subir à la matrice un shift, on remplace A par A−λ I. On peut prendre λ aléatoirement, ou bien mieux faire des itérations inverses sur A−λ I si λ est une estimation d’une valeur propre (voir les itérations inverses ci-dessous).

Lorsqu’on applique cette méthode a une matrice réelle, il peut arriver quíl y ait deux valeurs propres conjuguées de module maximal. On peut appliquer la méthode ci-dessus avec un shift complexe non réel, mais on doit alors travailler en arithmétique complexe ce qui est plus couteux. Le même type de raisonnement que ci-dessus montre que pour k grand, vk+2 est presque colinéaire à l’espace engendré par vk et vk+1, la recherche d’une relation avk+2+ b vk+1 + vk=0 permet alors de calculer les valeurs propres qui sont les deux racines de ax2+bx+1=0.

La convergence est de type série géométrique, on gagne le même nombre de décimales à chaque itération.

Applications :

Ceci peut par exemple servir à déterminer pour un polynôme P donné squarefree (de degré n et coefficient dominant pn) l’écart minimal entre 2 racines, on calcule
R:=normal(resultant(P,subst(P,x=x+y),x)/x^degree(P))
c’est un polynôme bicarré dont on cherche la plus petite racine en calculant le carré de la plus grande racine en module de numer(subst(R,y=1/sqrt(x))).

On peut obtenir un minorant à priori de cette plus petit racine en calculant

resultant(P,P′)=  ± pn2n−1 
 
1 ≤ i < j ≤ n
 (rirj)2 

on isole l’écart minimal au carré, on majore les autres carrés en majorant les racines, et on peut minorer le résultant à priori par 1 si P est à coefficients entiers.

7.2  Itérations inverses

La méthode précédente permet de calculer la valeur propre de module maximal d’une matrice. Pour trouver une valeur propre proche d’une quantité donnée x, on peut appliquer la méthode précédente à la matrice (AxI)−1 (en pratique on effectue LU sur AxI et on résoud (AxI)un+1=un). En effet, les valeurs propres de cette matrice sont les (xix)−1 dont la norme est maximale lorsqu’on se rapproche de xi. Attention à ne pas prendre x trop proche d’une valeur propre, car le calcul de (AxI)un+1=un est alors peu précis (la matrice étant mal conditionnée).

7.3  Elimination des valeurs propres trouvées

Si la matrice A est symétrique, et si en est un vecteur propre normé écrit en colonne, on peut considérer la matrice B=Axn en ent qui possède les mêmes valeurs propres et mêmes vecteurs propres que A avec même multiplicité, sauf xn qui est remplacé par 0. En effet les espaces propres de A sont orthogonaux entre eux, donc

Ben=xnen −xn en ent en = 0, Bek = xk ek − xn en ent ek = xk ek

On peut donc calculer la 2ème valeur propre (en valeur absolue), l’éliminer et ainsi de suite.

Si la matrice A n’est pas symétrique, il faut considérer B=Axn en fnt /en.fnfn est vecteur propre de At associé à xn. En effet fk.ej=0 si ij car Aek.fj=xkek.fj=ek.Atfj=xjekfj et donc fk.ek ≠ 0 (sinon ek est dans l’orthogonal de ℝn=Vect(f1,...,fn).

7.4  Décomposition de Schur

Il s’agit d’une factorisation de matrice sous la forme

A = P S P−1

P est unitaire et S diagonale supérieure. Existence (théorique) : on prend une valeur propre et un vecteur propre correspondant, puis on projette sur l’orthogonal de ce vecteur propre et on s’y restreint, on prend à nouveau une valeur propre et un vecteur propre correspondant, etc.

On peut approcher cette factorisation par un algorithme itératif qui utilise la factorisation QR d’une matrice quelconque comme produit d’une matrice unitaire par une matrice triangulaire supérieure à coefficients positifs sur la diagonale. On fait l’hypothèse que les valeurs propres de S sur la diagonale sont classées par ordre de module strictement décroissant |λ1|>|λ2|>...>|λn| (développement inspiré par Peter J. Olver dans le cas symétrique http://www.math.umn.edu/~olver/aims_/qr.pdf). On peut toujours s’y ramener quitte à remplacer A par A− α I. Posons A1=A, et par récurrence An=QnRn (avec Qn unitaire et R triangulaire supérieure à coefficients diagonaux positifs), An+1=RnQn. On a alors

 Ak=(Q1 R1) (Q1 R1) (Q1 R1) ... (Q1 R1) (Q1 R1
 =Q1 (R1 Q1) (R1 Q1) (R1 ... Q1) (R1 Q1R1 
 =Q1 (Q2 R2) (Q2 R2) .. (Q2 R2R1 
 =Q1 Q2 (R2 Q2R2 .. Q2 R2 R1 
 =Q1 Q2 (Q3 R3) ... Q3 R3 R2 R1 
 =Q1 ... Qk Rk ...R1

D’autre part A=PSP−1 donc Ak = P Sk P−1. Soit D la forme diagonale de S et U la matrice de passage S=UDU−1, où U est triangulaire supérieure et où on choisit la normalisation des coefficients sur la diagonale de U valant 1. On a donc

Ak = P U Dk U−1 P−1 

Ensuite, on suppose qu’on peut factoriser U−1P−1=LŨ sans permutations, donc qu’on ne rencontre pas de pivot nul, et quitte à multiplier les vecteurs unitaires de P−1 par une constante complexe de module 1 on peut supposer que les pivots sont positifs donc que Ũ a des coefficients positifs sur la diagonale, on a alors

Ak = P U Dk L Ũ = Q1 ... Qk Rk ...R1 

puis en multipliant par Ũ−1 |D|k

P U Dk L |D|k =  Q1 ... Qk Rk ...R1 Ũ−1 |D|k 

Rk ...R1 Ũ−1 |D|k est triangulaire supérieure à coefficients positifs sur la diagonale et Q1 ... Qk est unitaire. On regarde ensuite les entrées de la matrice Dk L |D|k, sous la diagonale elles convergent (géométriquement) vers 0, donc UDk L |D|k tend vers une matrice triangulaire supérieure dont les coefficients diagonaux valent ei karg(λj). On montre que cela entraine que Q1 ... Qk est équivalent à P(D/|D|)k

Q1 ... Qk ≈ P (D/|D|)k,      Rk ...R1 Ũ−1 |D|k ≈ (D/|D|)k UDk L |D|k 

Donc, Qk tend à devenir diagonale, et Rk Qk=Ak+1 triangulaire supérieure. De plus

A=Q1 A2 Q1−1 = ... = Q1 ... Qk Ak+1 (Q1 ... Qk)−1 

la matrice Ak+1 est donc semblable à A.

En pratique, on n’impose pas la positivité des coefficients diagonaux de R dans la factorisation QR, ce qui ne change évidemment pas le fait que Qk s’approche d’une matrice diagonale et Ak d’une matrice triangulaire supérieure (avec convergence à vitesse géométrique). On commence aussi par mettre la matrice A sous forme de Hessenberg (par conjugaison par des matrices de Householder), c’est-à-dire presque triangulaire supérieure (on autorise des coefficients non nuls dans la partie inférieure seulement sur la sous-diagonale, aij=0 si i>j+1). Cela réduit considérablement le temps de calcul de la décomposition QR, le produit RQ ayant encore cette propriété, une itération se fait en temps O(n2) au lieu de O(n3). Le calcul de RQ à partir de A est d’ailleurs fait directement, on parle d’itération QR implicite.

On utilise aussi des “shifts” pour accélerer la convergence, c’est-à-dire qu’au lieu de faire QR et RQ sur la matrice Ak on le fait sur Ak − αk I où λk est choisi pour accélerer la convergence vers 0 du coefficient d’indice ligne n colonne n−1 (idéalement il faut prendre αk proche de λn la valeur propre de module minimal, afin de minimiser |λn−αk|/|λn−1−αk|). En effet, si Ak − λk I = Qk Rk et Ak+1=RkQkk I alors :

 (A−α1 I) ... (A−αk I)=Q1R1(Q1R1−(α2−α1)I) ...(Q1R1−(αk−α1)I)
 =Q1 (R1Q1 − (α2−α1)I)R1  (Q1R1−(α3−α1)I) ...(Q1R1−(αk−α1)I)
 =Q1 (A2 − α1 I −(α2−α1)IR1  Q1(R1Q1−(α3−α1)IR1...(Q1R1−(αk−α1)I)
 =Q1 (A2 − α2 I) (A2 − α3 I) ... (A2 − αkIR1
 =...
 =Q1 ... Qk Rk ...R1

On peut aussi éliminer la dernière ligne et la dernière colonne de la matrice pour accélerer les calculs dès que le coefficient en ligne n colonne n−1 est suffisamment petit.

On remarque que pour une matrice réelle si on choisit des shifts conjugués, alors Q1...Qk Rk...R1 est réel. Or si QR=QR et si R est inversible

Q
−1 Q = 
R
 R−1 

On a donc une matrice symétrique (car Q−1=Qt) et triangulaire supérieure. On en déduit que Q−1 Q=D est diagonale, donc Q=Q D. On peut donc rendre Q réelle en divisant chaque colonne par un eiθ, et rendre R réelle en conjuguant par la matrice D. Mais ce procédé de retour au réel après élimination de 2 valeurs propres complexes conjuguées d’une matrice réelle se heurte à un problème de conditionnement parce que le choix d’un shift intéressant pour la convergence va rendre la matrice R proche d’une matrice non inversible (les deux derniers coefficients diagonaux de R sont proches de 0). On a alors seulement

Q
−1 Q R = 
R
 

Si on décompose Q−1 Q, R, R par blocs n−2,n−2, n−2,2, 2,n−2 et 2,2, on a



QQ11QQ12 
QQ21QQ22 
 



R11R12 
0R22 
 

=


QQ11 R11QQ21R11 
QQ11 R12 + QQ12 R22QQ21 R12 + QQ22 R22 
 

 =






R
11
R
12 
0
R
22 
 





Donc on a QQ11 =R11 R11−1. Comme Q est unitaire, QQ=Q−1 Q=Qt   Q est symétrique, donc QQ11 est diagonale puisque symétrique et triangulaire supérieure. On peut donc ramener Q11 et R11 en des matrices réelles. L’algorihtme des itérations QR implicites traite de manière efficace le cas des couples de valeurs propres complexes conjuguées ou plus généralement de clusters de valeurs propres, c’est l’algorithme de Francis (aussi appelé bulge chasing en anglais, qu’on pourrait traduire par “à la poursuite du bourrelet”, cela vient de la forme que prend la matrice après application d’un shift, elle a des entrées non nulles en première colonne plus bas que la sous-diagonale qui forment un bourrelet non nul, l’annulation de ces entrées par des transformations de Householder déplace le bourrelet sur la colonne suivante).

Revenons à la localisation des valeurs propres (une autre approche consiste à rechercher un rectangle du plan complexe stable par itérée de la méthode de Newton). On suppose qu’on a maintenant une matrice unitaire P et une matrice triangulaire supérieure S (aux erreurs d’arrondi près) telles que

P−1 A P = S 

Que peut-on en déduire ? On va d’abord arrondir P en une matrice exacte à coefficients rationnels, dont les dénominateurs sont une puissance de 2 (en fait c’est exactement ce que donne l’écriture d’un flottant en base 2, une fois ramené tous les exposants à la même valeur). On a donc une matrice Pe presque unitaire exacte et telle que

Se = Pe−1 A Pe  

est semblable à A, et presque triangulaire supérieure. (comme Pe est presque unitaire, sa norme et la norme de son inverse sont proches de 1 donc Se est proche de S, les coefficients de Se sont de la même taille que les coefficients de A : le changement de base est bien conditionné et c’est la raison pour laquelle on a choisi d’effectuer des transformations unitaires).

Notons µ1, ..., µn les coefficients diagonaux de Se, soit ε un majorant de la norme des coefficients sous-diagonaux de Se, et soit δ un minorant de l’écart entre 2 µj distincts. On a donc Se=U+EU est triangulaire supérieure, E est triangulaire inférieure avec des 0 sous la diagonale et des coefficients de module majorés par ε. Si ε est suffisamment petit devant δ, on va montrer qu’on peut localiser les valeurs propres de Se (qui sont celles de A) au moyen des µj.

En effet, fixons j et soit C un cercle de centre µ=µj et de rayon α ≤ δ/2. Si A est une matrice diagonalisable, on sait que

nombre de valeurs propres ∈ C =  
1
2iπ
  trace 
 


C
 (AzI)−1 

En prenant A=Se, et en écrivant

(SezI)−1 = (UzI+E)−1 = ( I + (UzI)−1E)−1(UzI)−1  

on développe le second terme si la norme de (UzI)−1E est strictement inférieure à 1

(SezI)−1 = (UzI)−1 −(UzI)−1E(UzI)−1 + (UzI)−1E(UzI)−1E(UzI)−1 + ...

puis on calcule la trace

trace(SezI)−1 = 
 
j
 (µjz)−1 + η 

avec

|η| ≤ 2π α || (UzI)−1 || 
|| (UzI)−1E ||
1−|| (UzI)−1E ||

Au final, le nombre de valeurs propres dans C est donné par

1 + η,    |η| ≤  α maxz ∈ C || (UzI)−1 || 
|| (UzI)−1E ||
1−|| (UzI)−1E ||
 

Il suffit donc que le max soit plus petit que 1 pour avoir l’existence d’une valeur propre et une seule de Se dans le cercle C (à distance au plus α de µ). Ce sera le cas si

ε ≤ 
1
2
 


δ
2 ||Se||
 


n−1



 
 
α
n−1

on choisit donc α pour réaliser l’égalité ci-dessus, sous réserve que δ ne soit pas trop petit, rappelons que α doit être plus petit ou égal à δ/2. Si δ est petit, il peut être nécessaire d’utiliser une précision plus grande pour les calculs de la décomposition de Schur en arithmétique flottante.

Typiquement, on peut espérer (pour un écart δ pas trop petit) pouvoir localiser les racines d’un polynôme de degré n par cette méthode avec précision b bits en O(n3 b2 + n2 b3) opérations pour le calcul de la décomposition de Schur en flottant (n3b2 pour Hessenberg initial puis n2b2 par itération et un nombre d’itérations proportionnel à b). Pour le calcul exact de Se, il faut inverser une matrice de taille n avec des coefficients de taille proportionnelle à b donc O(n4bln(n)) opérations (en modulaire, la taille des coefficients de l’inverse est O(nb ln(n))) puis calculer un produit avec une matrice n,n de coefficients de taille proportionnelle à b, soit O(n4b2ln(nb)) opérations. Asymptotiquement, on peut faire mieux avec des méthodes de multiplication et d’opérations matricielles par blocs. Pour éviter la perte d’un facteur n, on peut aussi ne pas faire de calculs en mode exact et controler les erreurs sur la matrice S. On peut regrouper les valeurs propres par “clusters” si elles sont trop proches à la précision de b bits. Pour la recherche des racines d’un polynôme P, on peut montrer, en calculant le résultant de P et de P′ qui est en module plus grand ou égal à 1, et en l’écrivant comme produit des carrés de différences des racines, et en majorant toutes les différences de racine sauf une à l’aide de la norme infinie de P, qu’il faut au pire b=O(n) bits pour séparer les racines).

8  Equations différentielles (résolution numérique)

8.1  Méthodes à un pas

On considère l’équation différentielle

y′=f(t,y),    t ∈ ℝ,    y(t) ∈ ℝdy(0)=y0

y(t) est la fonction inconnue cherchée et où f est une fonction régulière de t et y (par exemple C1 sur un domaine pour avoir existence et non recoupement des courbes intégrales dans ce domaine). On cherche à approcher numériquement y(t) pour t>0. On présente ici des méthodes de résolution numérique à un pas, dont le principe consiste à discrétiser l’intervalle [0,t] en des subdivisions en temps de petite taille [0,t1], [t1,t2], ..., [tn−1,tn=t]. Si yi est une valeur approchée de y(ti) la méthode à un pas se traduit par une relation de récurrence entre yi et yi+1 qui reflète une méthode d’intégration approchée de

y(ti+1)=y(ti)+
ti+1


ti
 f(t,y(t))  dt

Par exemple, la méthode d’Euler explicite utilise la méthode des rectangles à gauche

yi+1 = yi + (ti+1tif(ti,yi)=yi+hf(ti,yi)

h=ti+1ti (pour une méthode à pas variable, le pas h peut dépendre de i) alors que la méthode d’Euler implicite utilise la méthode des rectangles à droite

yi+1 = yi + (ti+1tif(ti+1,yi+1)=yi+hf(ti+h,yi+1)

cette dernière relation nécéssite de résoudre une équation pour déterminer yi+1 d’où son nom de méthode implicite. Plus généralement, la méthode de résolution revient à se donner une fonction Φ(t,y,h) et à poser :

yi+1=yi+hΦ(ti,yi,h)

pour la méthode d’Euler explicite, Φ(t,y,h)=f(t,y), pour la méthode d’Euler implicite, Φ s’obtient en résolvant une équation (par exemple avec la méthode du point fixe, pour h suffisamment petit).

Lorsqu’on compare la solution de l’équation et une valeur approchée obtenue par une méthode à un pas, il faut distinguer

Plus précisément, on a le résultat suivant :

Théorème 19   Soit y(t) la solution de y′=f(t,y), y(t0)=y0 sur [t0,T]. On considére une méthode de résolution à un pas :
yi+1=yi+hiΦ(ti,yi,hi)
Si la méthode est d’ordre p, i.e. si pour h=max(hi) l’erreur locale satisfait
|y(t) + hΦ(t,y(t),h)−y(t+h)|≤ Cp hp+1,  quad ∀ t ∈ [t0,T], h ≤ H 
et si la fonction Φ est lipschitzienne en y de constante Λ pour hH et y dans un voisinage de la solution y(t), i.e. si
|Φ(t,z,h)−Φ(t,y,h)| ≤ Λ |zy
alors l’erreur globale vérifie
|y(tn)−yn| ≤ hp 
Cp
Λ
(eΛ(tnt0)−1)

Par exemple, pour Euler explicite, Φ(t,y,h)=f(t,y), la constante Λ est la constante de Lipschitz de f, et on prendra pour C1 un majorant de 1/2|∂y f(t,y)| dans un voisinage de la solution y(t) pour t ∈ [t0,tn].

Pour prouver ce résultat, il faut déterminer comment se propagent les erreurs locales introduites à chaque pas. Par exemple, on a une erreur locale au pas 1 y(t1)−y1 donc une condition initiale modifiée pour le pas 2 y1 au lieu de y(t1). Cette erreur se propage au pas 2 en une erreur

y1y(t1)+h1(Φ(t1,y1,h1)−Φ(t1,y(t1),h1)|  ≤ (1+h1 Λ)|y1y(t1)| ≤ eh1 Λ |y1y(t1)|

De même aux pas suivants, donc au pas n l’erreur locale au pas 1 s’est propagée en une erreur inférieure ou égale à 

ehn−1 Λ... eh2 Λ eh1 Λ |y1y(t1)| = eΛ (tnt0) |y1y(t1)| ≤ Cph0p+1 eΛ (tnt0)

Il faut ensuite sommer les erreurs locales propagées de chaque pas

n−1
i=0
 Cp hip+1 eΛ (tnti)  ≤ Cp hp 
n−1
i=0
 hi eΛ (tnti) 

Comme eΛ(tnt) est positive décroissante sur [t0,tn], on peut majorer la somme par l’intégrale

Cp hp 
tn


t0
 eΛ (tnt)  dt 

d’où le résultat.

(Voir aussi Demailly ou Hairer)

8.2  Méthodes de Runge-Kutta (explicites)

Ce sont des méthodes explicites qui utilisent une méthode de Newton-Cotes pour approcher ∫f(t,y(t)) dt sur [ti,ti+1]. Pour simplifier les notations, notons ti=α, ti+1=β, on a alors

β


α
 f(t,y(t))  ≡ 
N
k=0
 ωk fk,yk)))

Pour estimer la valeur de fk,yk)), il est nécessaire d’approcher yk) ce qui se fait par une méthode de Newton-Cotes, en utilisant les estimations des valeurs des yj), j<k. On a donc des méthodes de Newton-Cotes avec un sous-ensemble croissant de points d’interpolation, donc pour chaque valeur de k une suite de coefficients ωj,k, j<k correspondant à la méthode de Newton-Cotes utilisée. Il faut aussi indiquer la valeur de αk en donnant un coefficient ck ∈ [0,1] tel que

αk = ti + ck (ti+1ti) = ti+ck h

En pratique on stocke un tableau dont les lignes donnent ck et les ckωj,k, j<k, et le calcul de yk) se fait ligne par ligne

yk) ≈ Yk=y0) + h 
k−1
j=0
 ckωj,k fj,yj)) 

. Par exemple pour la méthode d’Euler explicite, il y a deux lignes contenant 0 et un seul coefficient :

0: 
1:1

. Pour la méthode du point milieu, il y a trois lignes, la deuxière ligne exprime comment on estime y(ti+h/2), la troisième y(ti+1)=y(ti+h) :

0: 
  
1
2
:
1
2
 
1:01

on a donc

y(ti+
1
2
h) ≈ Y1=y(ti)+
1
2
hf(ti,y(ti))
y(ti+h) ≈ Y2=y(ti)+hf(ti+
h
2
,Y1)=y(ti)+hf(ti+
h
2
,y(ti)+
h
2
f(ti,y(ti)))

La suite des temps αk est croissante, mais pas forcément de manière stricte, on peut avoir αkk+1, la valeur de yk) n’étant pas estimée par la même méthode de Newton-Cotes que yk+1). La valeur des coefficients est ensuite déterminée pour obtenir un ordre le plus grand possible pour l’erreur locale (ce qui peut nécessiter la résolution de systèmes avec pas mal d’inconnues).

Ainsi, la méthode RK4 utilise le tableau suivant

0: 
1
2
:
1
2
 
1
2
:0
1
2
 
1:00
1:
1
6
1
3
1
3
1
6

Ce qui se traduit par

 Y1=
y(t0)+
h
2
f(t0,y0
Y2=
y(t0)+
h
2
f(t0+
h
2
,Y1
Y3=y(t0)+h f(t0+h,Y2
Y4=
y(t0) + 
h
6
 


f(t0,y(t0))+ 2f(t0+
h
2
,Y1)+2f(t0+
h
2
,Y2)+f(t0+h,Y3)


Les méthodes de Newton-Cotes utilisées sont les rectangles à gauche puis à droite pour estimer le point milieu, et la méthode de Simpson (en prenant la moyenne des deux estimations pour le point milieu). On peut montrer qu’elle est d’ordre 4 (erreur locale en O(h5))

Les méthodes de résolution numériques implémentées dans Xcas sont des méthodes explicites de Runge-Kutta emboitées avec pas adaptatif, (le pas adaptatif est calculé en estimant l’erreur avec 2 méthodes emboitées RK4 et Prince-Dormand, cf. Hairer).

9  Quelques références

Index

A  Développement de Taylor, séries entières, fonctions usuelles

Résumé: Séries entières. Calcul des fonctions transcendantes usuelles.

Soit f une fonction indéfiniment dérivable sur un intervalle I de ℝ et x0I. On peut alors effectuer le développement de Taylor de f en x0 à l’ordre n

Tn(f)(x)= f(x0) + (xx0f′(x0) + ... +  (xx0)n 
f[n](x0)
n!
 

et se demander si Tn(f) converge lorsque n tend vers l’infini, si la limite est égale à f(x) et si on peut facilement majorer la différence entre f(x) et Tn(f)(x). Si c’est le cas, on pourra utiliser Tn(f)(x) comme valeur approchée de f(x).

On peut parfois répondre à ces questions simultanément en regardant le développement de Taylor de f avec reste : il existe θ compris entre x0 et x tel que

Rn(x) := f(x)− Tn(f)(x) = (xx0)n+1
f[n+1](θ)
(n+1)!
 

C’est le cas pour la fonction exponentielle que nous allons détailler, ainsi que les fonctions sinus et cosinus.

A.1  La fonction exponentielle

Soit f(x)=exp(x) et x0=0, la dérivée n-ième de f est exp(x), donc Rn(x)=exp(θ)xn+1/(n+1)! avec θ compris entre 0 et x, ainsi si x est positif |Rn(x)| ≤ ex xn+1/(n+1)! et si x est négatif, |Rn(x)| ≤ xn+1/(n+1)!. Dans les deux cas, la limite de Rn est 0 lorsque n tend vers l’infini, car pour n ≥ 2x, on a

xn+1
(n+1)!
 = 
xn
n!
 
x
n+1
≤ 
1
2
xn
n!

on a donc pour tout x réel

ex = 
 
lim
n → +∞
 Tn(f)(x) = 
 
lim
n → +∞
 
n
k=0
 
xk
k!
k=0
xk
k!
 

Comment en déduire une valeur approchée de ex? Il suffira d’arrêter la sommation lorsque R:=xn+1/(n+1)! si x<0 ou lorsque R:=ex xn+1/(n+1)! si x>0 est inférieur à l’erreur absolue souhaitée, le plus tôt étant le mieux pour des raisons d’efficacité et pour éviter l’accumulation d’erreurs d’arrondi. Si on veut connaitre ex à une erreur relative ε donnée (par exemple ε=2−53 pour stocker le résultat dans un double) il suffit que R/ex < ε, donc si x est positif, il suffit que xn+1/(n+1)!<ε, on peut donc arrêter la sommation lorsque le terme suivant est plus petit que ε.

On observe que plus x est grand, plus n devra être grand pour réaliser le test d’arrêt, ce qui est facheux pour le temps de calcul. De plus, le résultat final peut être petit alors que les termes intermédiaires calculés dans la somme peuvent être grands, ce qui provoque une perte de précision relative, par exemple si on veut calculer e−10 ou plus généralement l’exponentielle d’un nombre négatif de grande valeur absolue.

Exercice : combien de termes faut-il calculer dans le développement de l’exponentielle de -10 pour que le reste soit plus petit que 2−53 ? Quel est la valeur du plus grand terme rencontré dans la suite ? Quelle est la perte de précision relative occasionné par cette méthode de calcul ?

On peut utiliser les propriétés de la fonction exponentielle pour éviter ce problème. Pour les nombres négatifs, on peut utiliser l’équation ex=1/ex (ne change pas l’erreur relative). Pour les grands réels, on peut utiliser e2x=(ex)2 (multiplie par 2 l’erreur relative). On peut aussi, si on connait une valeur approchée de ln(2), effectuer la division euclidienne de x par ln(2) avec reste symétrique :

x = a ln(2) + r,    a ∈ ℤ, |r| ≤ 
ln(2)
2
 

puis si r est positif, on somme la série de T(f)(r), si r est négatif, on calcule T(f)(−r) et on inverse, on applique alors :

ex = 2a er 

Il faut toutefois noter que ln(2) n’étant pas connu exactement, on commet une erreur d’arrondi absolu sur r d’ordre a η, où η est l’erreur relative sur ln(2), il faut donc ajouter une erreur d’arrondi relative de x/ln(2) η qui peut devenir grande si x est grand. Puis il faut ajouter la somme des erreurs d’arrondi due au calcul de er, que l’on peut minimiser en utilisant la méthode de Horner pour évaluer Tn(f)(r) (car elle commence par sommer les termes de plus haut degré qui sont justement les plus petits termes de la somme). Les coprocesseurs arithmétiques qui implémentent la fonction exponentielle ont un format de représentation interne des double avec une mantisse plus grande que celle des double (par exemple 64 bits au lieu de 53), et une table contenant des constantes dont ln(2) avec cette précision, le calcul de ex par cette méthode entraine donc seulement une erreur relative d’arrondi au plus proche sur le résultat converti en double (donc de 2−53).

Notons que en général x lui-même a déjà été arrondi ou n’est connu qu’avec une précision relative. Or si x>0 est connu avec une erreur relative de ε (donc une erreur absolue de ε |x|, alors

ex+ε |x|ex eε |x| 

donc on ne peut pas espérer mieux qu’une erreur relative de eε |x|−1 sur l’exponentielle de x. Si ε x est petit cette erreur relative (impossible à éviter, quel que soit l’algorithme utilisé pour calculer l’exponentielle) est d’ordre ε |x|. Si ε x est grand alors l’erreur relative devient de l’ordre de 1, et la valeur de l’exponentielle calculée peut être très éloignée de la valeur réelle! Notons que pour les double, il y aura dans ce cas débordement soit vers l’infini soit vers 0 (par exemple si x est supérieur à 709, l’exponentielle renvoie infini).

Exercice : refaire les mêmes calculs pour les fonction sinus ou cosinus. On utilise par exemple sin(x+π)=−sin(x), sin(−x)=−sin(x), sin(x)=cos(π/2−x) pour se ramener au calcul de sin(x) ou de cos(x) sur [0,π/4].

sin(x)=
n=0
(−1)n 
x2n+1
(2n+1)!
,    cos(x)=
n=0
(−1)n 
x2n
(2n)!
 

Cette méthode a toutefois ces limites, car il peut devenir impraticable de calculer la dérivée n-ième d’une fonction (par exemple avec tan(x)), et encore plus de la majorer. D’où l’intérêt de développer une théorie des fonctions qui sont égales à leur développement de Taylor à l’infini d’une part, et d’avoir d’autres méthodes pour majorer le reste, nous présentons ici le cas des séries alternées.

A.2  Séries entières.

Les séries de type prendre la limite lorsque n tend vers l’infini du développement de Taylor en x=0 sont de la forme

n=0
an xn := 
 
lim
 k → +∞
 
k
n=0
 an xnan=
f[n](0)
n!

On peut s’intéresser plus généralement à ∑n=0an xn lorsque an est un complexe quelconque, c’est ce qu’on appelle une série entière, on peut aussi les voir comme des polynômes généralisés.

S’il existe un point x0 tel que |an x0n| est borné (ce sera le cas en particulier si la série converge en x0), alors

|an xn| = |an x0n| |
x
x0
|n ≤  M |
x
x0
|n

la série converge donc en x si |x|<|x0| et on peut majorer le reste de la série au rang n par

|Rn| ≤ M 
 |
x
x0
|n+1
1−|
x
x0
|
 

la vitesse de convergence est donc du même type que pour le théorème du point fixe (le nombre de termes à calculer pour trouver une valeur approchée avec k décimales dépend linéairement k, les constantes sont d’autant plus grandes que |x| est grand).

Théorème 20   S’il existe un rang n0, un réel M>0 et un complexe x0 tels que pour n>n0, on ait :
|an x0|n ≤ M
alors la série converge pour |x|<|x0| et pour nn0, on a :
  |Rn| ≤ M 
 |
x
x0
|n+1
1−|
x
x0
|
      (14)

On en déduit qu’il existe un réel positif R≥ 0 éventuellement égal à +∞ tel que la série converge (la limite de la somme jusqu’à l’infini existe) lorsque |x|<R et n’existe pas lorsque |x|>R, ce réel est appelé rayon de convergence de la série. Par exemple ce rayon vaut +∞ pour l’exponentielle, le sinus ou le cosinus. Il est égal à 1 pour la série géométrique ∑xn (car elle diverge si |x|>1 et converge si |x|<1). On ne peut pas dire ce qui se passe génériquement lorsqu’on est à la limite, c’est-à-dire lorsque |x|=R (si R≠ +∞). Mais cela n’a en fait pas trop d’importance en pratique car même si la série converge, elle converge souvent trop lentement pour donner de bonnes approximations. En fait, la vitesse de convergence d’une série entière de rayon R≠ +∞ est en gros la même que celle d’une série géométrique de raison |x|/R.

Lorsque 2 séries ont un rayon de convergence non nul, alors on peut effectuer leur somme, leur produit comme des polynômes et la série somme/produit a un rayon de convergence au moins égal au plus petit des 2 rayons de convergence des arguments. On peut inverser une série entière non nulle en 0 en appliquant

(1+x)−1 = 1−x+x2x3+... 

et on obtient une série entière de rayon de convergence non nul. On peut aussi composer deux séries entières g et f en gf (avec les règles de calcul de composition des polynômes) si f(0)=0. On peut enfin dériver et intégrer une série entière terme à terme dans son rayon de convergence.

On dit qu’une fonction est développable en série entière en 0 si elle est égale à son développement de Taylor en 0 sommé jusqu’en l’infini dans un disque de centre 0 et de rayon non nul. Les fonctions exponentielle, sinus, cosinus sont donc développables en série entière en 0. La fonction tangente également car le dénominateur cosinus est non nul en 0, mais son rayon de convergence n’est pas l’infini et le calcul des an est assez complexe. La fonction (1+x)α est développable en séries entières pour tout α ∈ ℝ avec un rayon de convergence 1 (ou l’infini pour α entier positif).

(1+x)α= 1 + α x + 
α (α−1)
2!
 x2 + ... + 
α (α−1) ... (α −n +1)
n!
 xn + ...

Pour α=−1, c’est la série géométrique de raison −x, en effet si |x|<1 :

k
n=0
 (−x)n = 
1−(−x)k+1
1+x
  →k→ ∞ 
1
1+x

En intégrant par rapport à x, on obtient que ln(1+x) est développable en série entière en 0 de rayon de convergence 1 et

ln(1+x) = 
n=0
(−x)n+1
n+1
 

On peut calculer de manière analogue le développement en série entière de arctan(x) en iintégrant celui de 1/(1+x2), de même pour arccos(x) et arcsin(x) en intégrant celui de (1−x2)−1/2.

arctan(x)=
n=0
(−1)n 
x2n+1
2n+1
,

On peut donc calculer ln, arctan, ... par ces formules, mais il faut répondre à la question où arrête-t-on la somme pour obtenir une précision donnée? Dans le cas de ln(1+x), on pourrait répondre comme avec l’exponentielle en majorant la dérivée n+1-ième, mais ce n’est plus faisable pour arctan, arcsin, arccos. On va donner un autre critère qui ne nécessite pas de calculer cette dérivée mais utilise l’alternance des signes dans la somme.

A.3  Série alternée

Théorème 21   Soit Sn= ∑k=0n (−1)k uk la somme jusqu’au rang n d’une série de réels tels que la suite des uk décroit à partir d’un rang n0 et tend vers 0 lorsque k→ +∞. Alors Sn converge vers une limite S. Si nn0, la limite est comprise entre deux sommes partielles succesives Sn et Sn+1 et le reste est majoré par la valeur absolue du premier terme non sommé :
|Rn| ≤ |un+1|

Démonstration :
on montre que les suites vn=S2n et wn=S2n+1 sont adjacentes. On a

vn+1vnS2n+2S2n= (−1)2n+2 u2n+2 + (−1)2n+1 u2n+1 = u2n+2u2n+1 ≤ 0

donc vn est décroissante, de même wn est croissante, et vnwn=u2n+1 est positif et tend vers 0. On en déduit que vn et wn convergent vers la même limite S telle que vn>S>wn et les inégalités du théorème s’en déduisent.

Remarque
lorsqu’on utilise une suite alternée pour trouver une valeur approchée, il faut que un tende assez vite vers 0, sinon il y aura perte de précision sur la mantisse lorsqu’on effectuera u2nu2n+1. On sommera aussi les termes par ordre décroissant pour diminuer les erreurs d’arrondi.

A.4  La fonction logarithme

Si nous voulons calculer ln(1+x) pour x ∈ [0,1[ avec une précision ε, il suffit de calculer

n
k=0
 (−1)k 
xk+1
k+1

pour n tel que la valeur absolue du terme suivant soit plus petit que ε :

n  tel que  
xn+1
n+1
 < ε 

en effet, les signes sont alternés et la suite xk+1/k+1 décroit vers 0.

Si la suite décroit lentement vers 0, cette méthode est mauvaise numériquement et en temps de calcul car il y a presque compensation entre termes successifs donc perte de précision sur la mantisse et il y a beaucoup de termes à calculer. C’est le cas pour le logarithme, si x est voisin de 1, il faut calculer n termes pour avoir une précision en 1/n, par exemple 1 million de termes pour avoir une précision de 1e−6 (sans tenir compte des erreurs d’arrondi). Si x est proche de 1/2 il faut de l’ordre de −ln(ε)/ln(2) termes ce qui est mieux, mais encore relativement grand (par exemple 50 termes environ pour une précision en 1e−16, 13 termes pour 1e−4). On a donc intérêt à se ramener si possible à calculer la fonction en un x où la convergence est plus rapide (donc |x| le plus petit possible). Par exemple pour le calcul de ln(1+x) on peut :

Nous sommes donc en mesure de calculer précisément le logarithme ln(1+x) pour disons |x|<1/2. Pour calculer ln sur ℝ+, on se ramène à [1,2] en utilisant l’écriture mantisse-exposant, puis si x∈[3/2,2] on peut en prendre la racine carrée pour se retrouver dans l’intervalle souhaité. On peut aussi effectuer une division par √2.

Remarquons que si x est connu à une erreur relative ε près, comme

ln(x(1 ± ε))=ln(x) + ln(1 ± ε) 

ln(x) est connu à une erreur absolue de |ln(1 ± ε)| ≈ ε. Si ln(x) est proche de 0, on a une grande perte de précision relative.

Finalement, nous savons calculer ln et exp sous réserve d’avoir dans une table la valeur de ln(2). Pour calculer ln(2) précisément, on peut utiliser

ln(2)=−ln(1/2)=−ln(1−1/2) 

et le développement en série calculé en mode exact avec des fractions à un ordre suffisant, on majore le reste en utilisant que le terme général de la série ln(1+x) est borné par M=1 en x=1, donc d’après (14) :

|Rn| ≤ 
1
2n

(on peut même obtenir 1/(n2n) car on a besoin de M uniquement pour les termes d’ordre plus grand que n, on peut donc prendre M=1/n). Par exemple, pour avoir ln(2) avec une mantisse de 80 bits, on effectue une fois pour toutes avec un logiciel de calcul formel :
a:=sum((1/2)^k/k,k=1..80)|
puis la division en base 2 avec 81 bits de précision iquo(numer(a)*2^81,denom(a))

Exercice : pour les fonctions trigonométriques, il faut une méthode de calcul de π. On peut par exemple faire le calcul de 16 arctan(1/5)−4arctan(1/239) en utilisant le développement de la fonction arctan à un ordre suffisant.

A.5  Autres applications

On peut calculer certaines intégrales de la même manière, par exemple

1/2


0
 
1
1+x3

mais aussi des fonctions définies par des intégrales (cas de nombreuses fonctions spéciales).

A.5.1  Exemple : la fonction d’erreur (error fonction, erf)

Cette fonction est définie à une constante multiplicative près par :

f(x)=
x


0
 et2  dt 

On peut développer en séries entières l’intégrand (rayon de convergence +∞), puis intégrer terme à terme, on obtient

f(x)= 
+∞
n=0
 (−1)n 
x2n
n! (2n+1)

Ce développement converge très rapidement pour |x|≤ 1. Par contre, pour |x| grand, il faut calculer beaucoup de termes avant que le reste soit suffisamment petit pour être négligeable, et certains termes intermédiaires sont grands, ce qui provoque une perte de précision qui peut rendre le résultat calculé complètement faux. Contrairement à la fonction exponentielle, il n’y a pas de possibilité de réduire l’argument à une plage où la série converge vite. Il faut donc

Exercice : donner une valeur approchée de f(1) à 1e−16 près. Combien de termes faut-il calculer dans la somme pour trouver une valeur approchée de f(7) à 1e−16 près ? Comparer la valeur de f(7) et la valeur absolue du plus grand terme de la série, quelle est la perte de précision relative si on effectue les calculs en virgule flottante ? Combien de chiffres significatifs faut-il utiliser pour assurer une précision finale de 16 chiffres en base 10 ? Calculer le développement asymptotique en l’infini et déterminer un encadrement de f(7) par ce développement. Combien de termes faut-il calculer pour déterminer f(10) à 1e−16 près par le développement asymptotique et par le développement en séries ? Quelle est la meilleure méthode pour calculer f(10) ?

A.5.2  Recherche de solutions d’équations différentielles

On peut aussi appliquer les techniques ci-dessus pour calculer des solutions de certaines équations différentielles dont les solutions ne s’expriment pas à l’aide des fonctions usuelles, on remplace dans l’équation la fonction inconnue par son développement en séries et on cherche une relation de récurrence entre an+1 et an. Si on arrive à montrer par exemple qu’il y a une solution ayant un développement alternée, ou plus généralement, si on a une majoration |an+1/an|<C, alors le reste de la série entière est majoré par |anxn|/(1−|Cx|) lorsque |x|<1/C, on peut alors calculer des valeurs approchées de la fonction solution à la précision souhaitée en utilisant le développement en séries entières.

A.5.3  Exemple : fonctions de Bessel d’ordre entier

Soit m un entier positif fixé, on considère l’équation différentielle

x2 y′′ + x y′ + (x2m2)y=0 

dont on cherche une solution série entière y=∑k=0ak xk . En remplacant dans l’équation, si x est dans le rayon de convergence de la série (rayon supposé non nul), on obtient

k=0
k(k−1)ak xk + 
k=0
k ak xk  + 
k=0
(x2m2ak xk =0

soit encore

0=
k=0
(k2m2+x2ak xk  
 =
m2 a0 + (1−m2)a1 x + 
k=2
[(k2m2ak +ak−2]xk 

Par exemple, prenons le cas m=0. On a alors a0 quelconque, a1 nul et pour k≥ 2

ak = − 
ak−2
k2

Donc tous les a d’indice impair sont nuls. Les pairs sont non nuls si a0≠ 0, et ils sont de signe alterné. Soit x fixé, on observe que pour 2k > |x|,

|a2k x2k| < |a2k−2 x2k−2

donc la série ∑k=0ak xk est alternée à partir du rang partie entière de |x| plus un. Donc elle converge pour tout x (le rayon de convergence de y est +∞) et le reste de la somme jusqu’à l’ordre 2n est inférieur en valeur absolue à :

|R2n(x)| ≤ |a2n+2 x2n+2

Par exemple, pour avoir une valeur approchée à 1e−10 près de y(x) pour a0=1 et |x|≤ 1, on calcule y=∑k=02n ak xk , on s’arrête au rang n tel que

|a2n+2 x2n+2| ≤ |a2n+2| ≤ 10−10 

On remarque que :

a2n = 
(−1)n
22 42 ... (2n)2
 = 
(−1)n
22n n!2
 

donc n=7 convient.

Pour m ≠ 0, on peut faire un raisonnement analogue (les calculs sont un peu plus compliqués).

On a ainsi trouvé une solution y0 de l’équation différentielle de départ dont on peut facilement calculer une valeur approchée (aussi facilement que par exemple la fonction sinus pour |x| ≤ 1), on peut alors trouver toutes les solutions de l’équation différentielle (en posant y=y0 z et en cherchant z).

Exercice : faire de même pour les solutions de y′′−xy=0 (fonctions de Airy).

A.6  Développements asymptotiques et séries divergentes

Un développement asymptotique est une généralisation d’un développement de Taylor, par exemple lorsque le point de développement est en l’infini. De nombreuses fonctions ayant une limite en l’infini admettent un développement asymptotique en l’infini, mais ces développements sont souvent des séries qui semblent commencer par converger mais sont divergentes. Ce type de développement s’avère néanmoins très utile lorsqu’on n’a pas besoin d’une trop grande précision sur la valeur de la fonction.

Nous allons illustrer ce type de développement sur un exemple, la fonction exponentielle intégrale, définie à une constante près par

f(x)=
+∞


x
 
et
t
  dt 

On peut montrer que l’intégrale existe bien, car l’intégrand est positif et inférieur à et (qui admet −et comme primitive, cette primitive ayant une limite en +∞). Pour trouver le développement asymptotique de f en +∞, on effectue des intégrations par parties répétées, en intégrant l’exponentielle et en dérivant la fraction rationnelle

 f(x)=
[
et
t
]x+∞ − 
+∞


x
 
et
t2
  dt 
 =
ex
x
 − 
+∞


x
 
et
t2
  dt 
 =
ex
x
 − ([
et
t2
]x+∞ − 
+∞


x
 
−2et
t3
 =
ex
x
 − 
ex
x2
 + 
+∞


x
 
2et
t3
  dt 
 =... 
 =
ex


1
x
 − 
1
x2
 + 
2
x3
 + ... + 
(−1)n n!
xn+1



− 
+∞


x
 
(−1)n (n+1)!et
tn+2
  dt 
 =S(x) + R(x)

  S(x)=ex


1
x
 − 
1
x2
 + 
2
x3
 + ... + 
(−1)n n!
xn+1



,     R(x)=− 
+∞


x
 
(−1)n (n+1)!et
tn+2
  dt      (15)

Le développement en séries est divergent puisque pour x>0 fixé et n tendant vers l’infini

 
lim
n→ +∞
 
n!
xn+1
 = +∞

mais si x est grand, au début la série semble converger, de manière très rapide :

1
x
 >> 
1
x2
 >> 
2
x3
 

On peut utiliser S(x) comme valeur approchée de f(x) pour x grand si on sait majorer R(x) par un nombre suffisamment petit. On a

R(x) | ≤ 
+∞


x
 
(n+1)!et
xn+2
(n+1)!ex
xn+2
 

On retrouve une majoration du type de celle des séries alternées, l’erreur relative est inférieure à la valeur absolue du dernier terme sommé divisé par ex/x. Pour x fixé assez grand, il faut donc trouver un rang n, s’il en existe un, tel que (n+1)!/xn+1<є où є est la précision relative que l’on s’est fixée. Par exemple, si x≥ 100, n=11 convient pour є=12!/10012=5e−16 (à peu près la précision relative d’un “double”). Ceci permet d’avoir une approximation de la fonction avec une bonne précision et peu de calculs, mais contrairement aux séries entières, il n’est pas possible d’améliorer cette précision de manière arbitraire en poussant le développement plus loin, il y a une précision maximale possible (qui dépend de x).

Ce type de développement asymptotique peut être effectué pour d’autres fonctions du même type, par exemple

+∞


x
 et2  dt,    
+∞


x
 
sin(t)
t
  dt,    ... 

Digression: calcul approché de la constante d’Euler γ
On peut montrer que

 
 
lim
n→ +∞
 un,    un=
n
k=1
1
k
 − ln(n)      (16)

existe (par exemple en cherchant un équivalent de un+1un qui vaut −1/2n2) et on définit γ comme sa limite. Malheureusement, la convergence est très lente et cette définition n’est pas applicable pour obtenir la valeur de γ avec une très grande précision. Il y a un lien entre γ et la fonction exponentielle intégrale, plus précisément lorsque x→ 0, f(x) admet une singularité en −ln(x), plus précisément f(x)+ln(x) admet un développement en séries (de rayon de convergence +∞), car :

 f(x)+ln(x)=
1


x
et−1
t
  dt + 
+∞


1
 
et
t
  dt 
 =
1


0
et−1
t
  dt + 
+∞


1
 
et
t
  dt − 
x


0
 
et−1
t
  dt

Que vaut la constante du membre de droite :

C=
1


0
(et−1)
1
t
  dt + 
+∞


1
 et 
1
t
  dt 

Il se trouve que C=−γ (voir plus bas une démonstration condensée) et donc :

  γ= 
x


0
 
1−et
t
  dt −f(x)−ln(x)     (17)

Pour obtenir une valeur approchée de γ, il suffit donc de prendre un x assez grand pour pouvoir calculer f(x) par son développement asymptotique à la précision requise, puis de calculer l’intégrale du membre de droite par le développement en séries en x=0 (en utilisant une précision intermédiaire plus grande puisque ce développement en séries va sembler diverger au début avant de converger pour n suffisamment grand). Par exemple, on pose x=13, on calcule f(13) par (15) avec n=13 (qui correspond au moment où le terme général de la série est minimum puisque le rapport de deux termes successifs est en n/x) et une erreur absolue inférieure à e−13 13!/1314=4e−12

f(13) ≈ exp(-13)*sum((-1)^n*n!/13.^(n+1),n=0..13)

puis on remplace dans (17), avec

x


0
 
1−et
t
  dt = 
n=0
 (−1)n 
xn+1
(n+1) (n+1)!

dont on obtient une valeur approchée, en faisant la somme jusqu’au rang 49 (pour lequel le terme général est de l’ordre de 1e-12), le reste de cette somme R50 est positif et est inférieur à (-1)^50*13.^51/51/51!) qui est de l’ordre de 8e-12

evalf(sum((-1)^n*13^(n+1)/(n+1)/(n+1)!,n=0..49))

La somme argument de evalf étant exacte, il n’y a pas de problèmes de perte de précision, on peut aussi faire les calculs intermédiaires en arithmétique approchée, on doit alors prendre 4 chiffres significatifs de plus pour tenir compte de la valeur du plus grand terme sommé dans la série, terme que l’on détermine par exemple par

seq(13.^(n+1)/(n+1)/(n+1)!,n=0..20)

ce terme vaut 13^11/11/11! soit 4000 environ)

Digits:=16; sum((-1)^n*13.^(n+1)/(n+1)/(n+1)!,n=0..49)

On obtient finalement comme valeur approchée de γ

-exp(-13)*sum((-1)^n*n!/13.^(n+1),n=0..13)-ln(13)+
sum((-1)^n*13^(n+1)/(n+1)/(n+1)!,n=0..49)

soit 0.577215664897 avec une erreur inférieure à 1.2e-11. Bien entendu, cette méthode est surtout intéressante si on veut calculer un grand nombre de décimales de la constante d’Euler, sinon on peut par exemple appliquer la méthode d’accélération de Richardson à la suite convergente (16) qui définit γ ou d’autres méthodes d’accélération (en transformant par exemple la série en série alternée). On calcule alors de deux manières différentes f(x) pour x plus grand (déterminé par la précision qu’on peut obtenir par le développement aymptotique de f).

On peut calculer π de la même manière avec le développement en séries et asymptotique de la fonction sinus intégral (on remplace exponentielle par sinus dans la définition de f) et l’égalité (dont un schéma de preuve est aussi donné plus bas)

 
+∞


0
 
sin(t)
t
  dt = 
π
2
    (18)

Calcul de C (et preuve de (18)):
Pour cela on effectue une intégration par parties, cette fois en intégrant 1/t et en dérivant l’exponentielle (moins 1 dans la première intégrale).

 C=
1


0
(et−1)
1
t
  dt + 
+∞


1
 et 
1
t
  dt
 =
[(et−1)ln(t)]01 +
1


0
 ln(tet  dt + [et ln(t)]1+∞ +
+∞


1
 ln(tet  dt 
 =
+∞


0
 ln(tet  dt

Pour calculer cette intégrale, on utilise l’égalité (qui se démontre par récurrence en faisant une intégration par parties) :

n!= 
+∞


0
tn et  dt 

On va à nouveau intégrer par parties, on intègre un facteur multiplicatif 1 et on dérive l’intégrand, on simplifie, puis on intègre t et on dérive l’autre terme, puis t2/2, etc.

 C=
[tet ln(t)]0+∞ − 
+∞


0
 t et(
1
t
−ln(t))  dt 
 =
0 − 
+∞


0
 et  dt + 
+∞


0
 t et ln(t)  dt 
 =
−1 + [
t2
2
et ln(t)]0+∞  − 
+∞


0
 
t2
2
 et(
1
t
−ln(t))  dt 
 =
−1 − 
+∞


0
 
t
2
 et +  
+∞


0
 
t2
2
 et ln(t)  dt 
 =
−1 − 
1
2
 +  
+∞


0
 
t2
2
 et ln(t)  dt 
 =...
 =
−1 − 
1
2
 − ... −  
1
n
 + 
+∞


0
 
tn
n!
 et ln(t)  dt 
 =
−1 − 
1
2
 − ... −  
1
n
 + ln(n) + In

In=
+∞


0
 
tn
n!
 et (ln(t)−ln(n))  dt 

Pour déterminer In on fait le changement de variables t=nu

 In=
+∞


0
 
(nu)n
n!
 enu ln(un du 
 =
nn+1
n!
 
+∞


0
 en(ln(u)−u) ln(u)  du 

Or en faisant le même changement de variables t=nu :

n!= 
+∞


0
tn et  dt = nn+1 
+∞


0
 en(ln(u)−u)  du

Donc

In
+∞


0
 en(ln(u)−u) ln(u)  du
+∞


0
 en(ln(u)−u)  du
  

Lorsque n tend vers l’infini, on peut montrer que In → 0, en effet les intégrales sont équivalentes à leur valeur sur un petit intervalle autour de u=1, point où l’argument de l’exponentielle est maximal, et comme l’intégrand du numérateur a une amplitude ln(u) qui s’annule en u=1, il devient négligeable devant le dénominateur. Finalement on a bien C=−γ.

On peut remarquer qu’en faisant le même calcul que C mais en remplacant et par e−α t pour ℜ(α)>0, donne limIn=−ln(α) (car le point critique où la dérivée de la phase s’annule est alors 1/α). Ceci peut aussi se vérifier pour α réel en faisant le changement de variables α t=u

1


0
(e−α t−1)
1
t
  dt + 
+∞


1
 e−α t 
1
t
  dt  = −γ −ln(α) 

En faisant tendre α vers −i, −ln(α) tend vers ln(i)=iπ/2 et on obtient

1


0
(eit−1)
1
t
  dt + 
+∞


1
 ei t 
1
t
  dt  = −γ + i 
π
2
 

dont la partie imaginaire nous donne (18), et la partie réelle une autre identité sur γ faisant intervenir la fonction cosinus intégral.

B  La moyenne arithmético-géométrique.

B.1  Définition et convergence

Soient a et b deux réels positifs, on définit les 2 suites

  u0=av0=b,    un+1=
un+vn
2
vn+1=
unvn
      (19)

On va montrer que ces 2 suites sont adjacentes et convergent donc vers une limite commune notée M(a,b) et il se trouve que la convergence est très rapide, en raison de l’identité :

  un+1vn+1=
1
2
(
un
vn
)2 =
1
2(
un
+
vn
)2
(unvn)2     (20)

la convergence est quadratique.

On suppose dans la suite que ab sans changer la généralité puisque échanger a et b ne change pas la valeur de un et vn pour n>0. On a alors unvn (d’après (20) pour n>0) et un+1un car

un+1un=
1
2
(vnun) ≤ 0

et vn+1=√unvn ≥ √vnvn=vn. Donc (un) est décroissante minorée (par v0), (vn) est croissante majorée (par u0), ces 2 suites sont convergentes et comme un+1=un+vn/2, elles convergent vers la même limite l qui dépend de a et b et que l’on note M(a,b). On remarque aussi que M(a,b)=bM(a/b,1)=aM(1,b/a).

Précisons maintenant la vitesse de convergence lorsque ab>0. On va commencer par estimer le nombre d’itérations nécessaires pour que un et vn soient du même ordre de grandeur. Pour cela, on utilise la majoration

ln(un+1)−ln(vn+1) ≤ ln(un)−ln(vn+1) = 
1
2
(ln(un)−ln(vn)) 

donc

ln
un
vn
 = ln(un)−ln(vn) ≤ 
1
2n
 (ln(a)−ln(b)) = 
1
2n
 ln
a
b
 

Donc si n ≥ ln( ln(a/b)/m)/ln(2) alors lnun/vnm (par exemple, on peut prendre m=0.1 pour avoir un/vn ∈ [1,e0.1]). Le nombre minimum d’itérations n0 est proportionnel au log du log du rapport a/b. Ensuite on est ramené à étudier la convergence de la suite arithmético-géométrique de premiers termes a=un0 et b=vn0 et même en tenant compte de M(a,b)=aM(1,b/a) à a=1 et b=vn/un donc 0≤ ab ≤ 1−e−0.1. Alors l’équation (20) entraine

un+1vn+1 ≤ 
1
8
(unvn)2 

puis (par récurrence)

0 ≤ unvn ≤ 
1
82n−1
(ab)2n 

Donc comme M(a,b) est compris entre vn et un, l’erreur relative sur la limite commune est inférieure à une précision donnée є au bout d’un nombre d’itérations proportionnel au ln(ln(1/є)).

Typiquement dans la suite, on souhaitera calculer M(1,b) avec b de l’ordre de 2n en déterminant n chiffres significatifs, il faudra alors O(ln(n)) itérations pour se ramener à M(1,b) avec b∈ [e−0.1,1] puis O(ln(n)) itérations pour avoir la limite avec n chiffres significatifs.

Le cas complexe
On suppose maintenant que a, b ∈ ℂ avec ℜ(a)>0, ℜ(b)>0. On va voir que la suite arithmético-géométrique converge encore.
Étude de l’argument
On voit aisément (par récurrence) que ℜ(un)>0 ; de plus ℜ(vn) > 0 car par définition de la racine carrée ℜ(vn)≥ 0 et est de plus non nul car le produit de deux complexes d’arguments dans ]−π/2,π/2[ ne peut pas être un réel négatif. On en déduit que arg(un+1)=arg(un+vn) se trouve dans l’intervalle de bornes arg(un) et arg(vn) et que arg(vn+1)=1/2(arg(un)+arg(vn)) donc

| arg(un+1−arg(vn+1) | ≤ 
1
2
|arg(un)−arg(vn)| 

Après n itérations, on a

|arg(un)−arg(vn)| ≤ 
π
2n
 

Après quelques itérations, un et vn seront donc presque alignés. Faisons 4 itérations. On peut factoriser par exemple vn et on est ramené à l’étude de la suite de termes initiaux a=un/vn d’argument arg(un)−arg(vn) petit (inférieur en valeur absolue à π/16) et b=1. On suppose donc dans la suite que

|arg(
un
vn
)| ≤ 
π/16
2n
 

Étude du module
On a :

un+1
vn+1
1
2








un
vn
+
1
un
vn








Posons un/vnn eiθn, on a :

|
un+1
vn+1
|
=
1
2





ρn
 eiθn/2
1
ρn
 eiθn/2 




 =
1
2
 




(
ρn
1
ρn
)cos
θn
2
i (
ρn
− 
1
ρn
)sin
θn
2
 




 =
1
2
 
 (
ρn
1
ρn
)2cos2
θn
2
+ (
ρn
− 
1
ρn
)2sin2
θn
2
 
 
 =
1
2
 
 ρn
1
ρn
 +2cosθn 

Si ρ désigne le max de ρn et 1/ρn, on a alors la majoration

|
un+1
vn+1
| ≤ 
1
2
 
 ρ + ρ + 2 ρ 
  = 
ρ
 

donc en prenant les logarithmes

  lnρn+1 ≤ 
1
2
  lnρ=
1
2
  |lnρn|      (21)

On rappelle qu’on a la majoration

|arg(
un
vn
)| = |θn| ≤ 
π/16
2n
 ≤ 
1
2n+1
 

qui va nous donner la minoration de ρn+1

ρn+1=|
un+1
vn+1
|
=
1
2
 ρn
1
ρn
 +2 − 2 (1−cosθn
 
 =
1
2
 
 ρn
1
ρn
 +2 − 4 sin2 (
θn
2
 
 
1
2
 
 ρn
1
ρn
 +2 − θn2
 
 
1
2
 
 ρn
1
ρn
 +2
 × 
1 − 
θn2
ρn
1
ρn
 +2
 
 
1
2
 
 
1
ρ
 + 
1
ρ
 +2
1
ρ
 × 
1 − 
θn2
4
 
 
1
ρ
 
1 − 
θn2
4
 
 
1
ρ
 
1 − 
1
4 × 22n+2

en prenant les log et en minorant ln(1−x) par −2x

lnρn+1 ≥ 
1
2
 (−|lnρn|+ln(1 −
1
4 × 22n+2
 )) ≥ −
1
2
 (|lnρn|+
1
22n+3
 )  

Finalement avec (21)

|lnρn+1| ≤ 
1
2
 (|lnρn|+
1
22n+3
 ) 

On en déduit

|lnρn| ≤ 
1
2n
 lnρ0 + 
1
2n+3
 + ... +
1
22n+1
 +  
1
22n+2
1
2n
 lnρ0 + 
1
2n+2
 

La convergence du ln(un/vn) vers 0 est donc géométrique, donc un et vn convergent quadratiquement.

B.2  Lien avec les intégrales elliptiques

Le calcul de la limite commune des suites un et vn en fonction de a et b n’est pas trivial au premier abord. Il est relié aux intégrales elliptiques, plus précisément on peut construire une intégrale dépendant de deux paramètres a et b et qui est invariante par la transformation un,vnun+1,vn+1 (19)

I(a,b)=
+∞


−∞
  
dt
(a2+t2)(b2+t2)

On a en effet

I(
a+b
2
,
ab
) = 
+∞


−∞
  
du
((
a+b
2
)2+u2)(ab+u2)
 

On pose alors

u=
1
2
 (t
ab
t
),    t>0 

tu est une bijection croissante de t∈]0,+∞[ vers u ∈ ]−∞,+∞[, donc

I(
a+b
2
,
ab
)
=
+∞


0
  
dt/2(1+ab/t2)
((
a+b
2
)2+1/4(tab/t)2)(ab+1/4(tab/t)2)
 =
+∞


0
 
dt
(a2+t2)(b2+t2)
I(a,b)

On note au passage que I est définie si a,b ∈ ℂ vérifient ℜ(a)>0, ℜ(b)>0, on peut montrer que la relation ci-dessus s’étend (par holomorphie).

Lorsque a=b=l (par exemple lorsqu’on est à la limite), le calcul de I(l,l) est explicite

I(l,l)=
+∞


−∞
  
dt
(l2+t2)
 = 
π
l

donc

I(a,b)=I(M(a,b),M(a,b))=
π
M(a,b)

On peut transformer I(a,b) en posant t=bu

I(a,b)=2
+∞


0
  
du
(a2+b2u2)(1+u2)
2
a
 
+∞


0
  
du
(1+(b/a)2u2)(1+u2)
 

Puis en posant u=tan(x) (du=(1+u2) dx)

I(a,b)=
2
a
 
π
2


0
 
1+tan(x)2
1+(b/a)2tan(x)2
  dx 

et enfin en posant tan2(x)=sin(x)2/1−sin(x)2

I(a,b)= 
2
a
 
π
2


0
  
 
1
1−(1−
b2
a2
)sin(x)2
 
  dx

Si on définit pour m<1

K(m)=
π
2


0
 
dx
1−m sin(x)2
 

alors on peut calculer K en fonction de I, en posant m=1−b2/a2 soit b2/a2=1−m

K(m)=
a
2
 I(a,a
1−m
)=
a
2
π
M(a,a
1−m
)
=
π
2M(1,
1−m
)
 

d’où l’on déduit la valeur de l’intégrale elliptique en fonction de la moyenne arithmético-géométrique :

  K(m)=
π
2


0
 
dx
1−m sin(x)2
π
2M(1,
1−m
)
      (22)

Dans l’autre sens, pour x et y positifs

K( (
xy
x+y
)2 )=  
π
2M(1,
1−(
xy
x+y
)2
)
=  
π
2M(1,
2
x+y
xy
)
π
2
x+y
 M(
x+y
2
,
xy
π
4
 
x+y
M(x,y)

et finalement

M(x,y)=
π
4
 
x+y
 K


(
xy
x+y



2



 
 )

B.3  Application : calcul efficace du logarithme.

On peut utiliser la moyenne arithmético-géométrique pour calculer le logarithme efficacement, pour cela on cherche le développement asymptotique de K(m) lorsque m tend vers 1. Plus précisément, on va poser 1−m=k2 avec k ∈ ]0,1], donc

K(m)= 
π
2


0
 
dx
1−(1−k2) sin(x)2
=
π
2


0
 
dy
1−(1−k2) cos(y)2

en posant y=π/2−x, et

K(m)=
π
2


0
 
dy
sin(y)2+k2 cos(y)2

la singularité de l’intégrale pour k proche de 0 apparait lorsque y est proche de 0. Si on effectue un développement de Taylor en y=0, on trouve

sin(y)2+k2 cos(y)2 = k2 + (1−k2y2 + O(y4)

Il est donc naturel de comparer K(m) à l’intégrale

J=
π
2


0
 
dy
k2 + (1−k2y2
 

qui se calcule en faisant par exemple le changement de variables

y=
k
1−k2
 sinh(t)

ou directement avec Xcas,

supposons(k>0 && k<1);
J:=int(1/sqrt(k^2+(1-k^2)*y^2),y,0,pi/2)

qui donne après réécriture :

  J
1
1−k2





ln


π
k



+ ln




1
2
 




 1−k2 +4 
k2
π2
  +
1−k2
 














    (23)

et on peut calculer le développement asymptotique de J en 0

series(J,k=0,5,1)

qui renvoie :

J =ln


π
k



+O


−1
ln(k)



5



 
)

on peut alors préciser ce développement par

series(J+ln(k)-ln(pi),k=0,5,1)

qui renvoie (après simplifications et où la notation Õ peut contenir des logarithmes)




1
π2
 + 
ln(π)−ln(k)−1
2



k2 +  Õ(k4)

donc

J=−ln(k)+ln(π)+


1
π2
 + 
ln(π)−ln(k)−1
2



k2 + Õ(k4)     (24)

Examinons maintenant KJ, il n’a plus de singularité en y=0, et il admet une limite lorsque k → 0, obtenue en remplacant k par 0

(KJ)|k=0 = 
π
2


0
 


1
sin(y)
1
y



 dy = 


ln


tan


y
2






− ln(y


π
2



0
 = ln(
4
π
)

D’où pour K

Kk → 0 = ln


4
k



O


−1
ln(k)



5



 
)

Pour préciser la partie du développement de K en puissances de k, nous allons majorer KJ−ln(4/π), puis J−ln(π/k). Posons

A=sin(y)2+k2 cos(y)2,    B=y2+(1−y2)k2

Majoration de KJ−ln(4/π)
L’intégrand de la différence KJ−ln(4/π) est

     
1
A
 − 
1
B
 − 


1
sin(y)
1
y
 


=
B
A
A
 
B
 −
y−sin(y)
ysin(y)
 
    (25)
 =
BA
A
 
B
 (
A
+
B
)
 −
y−sin(y)
ysin(y)
  
    (26)
 =
 
(y2−sin(y)2)(1−k2)
A
 
B
 (
A
+
B
)
− 
y−sin(y)
ysin(y)
  
    (27)

Soit

  KJ−ln(
4
π
)= 
π
2


0
 
(y−sin(y))[(1−k2)y sin(y)(y+sin(y))−
AB
(
A
+
B
)]
A
 
B
 (
A
+
B
)ysin(y)
    (28)

On décompose l’intégrale en 2 parties [0,k] et [k,π/2]. Sur [0,k] on utilise (26), on majore chaque terme séparément et on minore A et B par

A=k2+(1−k2)sin(y)2 ≥ k2,    B=k2+(1−k2)y2 ≥ k2

Donc

k


0
 |
k


0
 
|BA|
2k3
  dy + 
k


0
 ( 
1
sin(y)
1
y
 )   dy 
 
k


0
 
y2−sin(y)2
2k3
  dy + ln(tan(
k
2
)) −ln(
k
2
 
1
3
 k3+
−1
2
 k+
1
4
 sin(2 k)
k3
  + ln(sin(
k
2
)) −ln(
k
2
) − ln(cos(
k
2
))
 
1
3
 k3+
−1
2
 k+
1
4
 (2k
8k3
6
+
32k5
5!
k3
 − ln(cos(
k
2
)) 
 
k2
30
−  ln(1− 
1
2!



k
2



2



 
 
k2
30
 +
k2
4

Sur [k,π/2], on utilise (28) et on minore A et B par

A=sin(y)2+k2 cos(y)2 ≥ sin(y)2,    B=y2+(1−y2)k2 ≥ y2

on obtient

π
2


k
 | ≤  
π
2


k
(y−sin(y))|C|</