Mat 406

Bernard.Parisse@univ-grenoble-alpes.fr

2022

N.B. : la version HTML de ce document est interactive, vous pouvez exécuter certains champs d’entrée donnés en exemple en les modifiant éventuellement, il suffit de cliquer sur le bouton exe. Il est recommandé d’utiliser Firefox ou un navigateur compatible supportant MathML.

Table des matières

1  Présentation du module

Dans ce module, on introduira moins de notions nouvelles que dans d’autres modules de mathématiques, par contre on insistera sur le calcul effectif, si possible efficace, et sur le controle de la précision des résultats, ceci explique la part importante consacrée aux TP (18h de cours en 12 séances, 18h de TD en 12 séances et 19.5h de TP en 13 séances). On présentera par exemple des méthodes de calcul des fonctions usuelles (racine carrée, trigonométriques, ...), il s’agira non seulement de savoir calculer une valeur numérique, mais aussi de pouvoir majorer l’écart entre la valeur trouvée et la valeur exacte, en utilisant des théorèmes du cours. Les calculs se feront dans la mesure du possible sur ordinateur ou sur calculatrices.

Les thèmes abordés seront :

  1. Calcul exact et approché, représentation des données
  2. Suites récurrentes, méthode du point fixe, de Newton.
  3. Séries de Taylor et approximation des fonctions usuelles
  4. Arithmétique des polynômes (PGCD, Bézout, factorisation, décomposition en éléments simples)
  5. Interpolation polynomiale
  6. Intégration approchée.
  7. Algèbre linéaire.

L’ordre des deux premiers thèmes sera inversé en cours.

L’évaluation se fait sur :

Les calculatrices sont autorisées au DS et à l’examen final (prêt possible pour le semestre).

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

Résumé:
Types de base : entier machine, entier long, flottant machine et multiprécision (Base 2, base 10, BCD).
Types composés : complexes, polynomes (représentation dense/creuse), symboles, listes (vecteurs, matrices), expressions, fonctions.
Erreur relative, erreur absolue, erreur d’arrondi, +/-, */%
Algorithme, complexité, exemple puissance modulaire, algorithme de Horner.

Les principaux ensembles de nombres en mathématiques sont les entiers positifs \mathbb{N} et relatifs \mathbb{Z}, les rationnels \mathbb{Q}, les réels \mathbb{R} et les complexes \mathbb{C}. 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 aa et bb sont deux entiers, a0,b>0a \geq 0, b>0, il existe un unique couple (q,r)(q,r) tel que a=bq+r,r[0,b[ a = bq +r , \quad r \in [0, b[

Preuve : On prend pour qq le plus grand entier tel que abq0a-bq \geq 0.
Exemple :

La division euclidienne permet d’écrire un nombre entier, en utilisant une base bb et des caractères pour représenter les entiers entre 0 et b1b-1. Nous écrivons les nombres entiers en base b=10b=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 bb (par exemple b=2b=2), on effectue des divisions euclidienne successives par bb 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 dd à partir de son écriture d n...d 0d_n...d_0, on traduit les divisions euclidiennes successives en d = (...((d nb+d n1)b+d n2)...+d 1)b+d 0 = d nb n+d n1b n1+...+d 0 \begin{matrix} d &=&( ... ((d_n b +d_{n-1})b + d_{n-2})...+d_1)b+d_0\\ &=& d_n b^n + d_{n-1} b^{n-1} + ... + d_0 \end{matrix} 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=2 4+2 3+125=2^4+2^3+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 aa et bb 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 [2 31,2 311][-2^{31},2^{31}-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, 2 30+2 30+2 30+2 302^{30}+2^{30}+2^{30}+2^{30} 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 410^{-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 mm et l’exposant ee. Si on écrit les nombres en base bb, la mantisse mm s’écrira avec un nombre nn fixé de chiffres (ou de bits en base 2), donc m[0,b n[m \in [0,b^n[. Soit un réel xx représenté par x=mb e,m[0,b n[ x=mb^e, \quad m \in [0,b^n[ Si m[0,b n1[m\in [0,b^{n-1}[, alors on peut aussi écrire x=mb e1x=m' b^{e-1} avec m=mb[0,b n[m'=mb \in [0,b^n[, quelle écriture faut-il choisir? Intuitivement, on sent qu’il vaut mieux prendre mm' le plus grand possible, car cela augmente le nombre de chiffres significatifs (alors que des 0 au début de mm ne sont pas significatifs). Ceci est confirmé par le calcul de l’erreur d’arrondi pour représenter un réel. En effet, si xx est un réel non nul, il ne s’écrit pas forcément sous la forme mb emb^e, on doit l’arrondir, par exemple au plus proche réel de la forme mb emb^e. La distance de xx à ce réel est inférieure ou égale à la moitié de la distance entre deux flottants consécutifs, mb emb^e et (m+1)b e(m+1)b^e, donc l’erreur d’arrondi est inférieure ou égale à b e/2b^e/2. Si on divise par xmb ex \geq mb^e, on obtient une erreur relative d’arrondi majorée par 1/(2m)1/(2m). On a donc intérêt à prendre mm le plus grand possible pour minimiser cette erreur. Quitte à mulitplier par bb, on peut toujours se ramener (sauf exceptions, cf. ci-dessous), à m[b n1,b n[m \in [b^{n-1},b^n[, on a alors une erreur d’arrondi relative majorée par 12b n1 \frac{1}{2b^{n-1}}

On appelle flottant normalisé un flottant tel que m[b n1,b n[m \in [b^{n-1},b^n[. Pour écrire un réel sous forme de flottant normalisé, on écrit le réel en base bb, et on déplace la virgule pour avoir exactement nn 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 532^{-53} (b=2b=2 et n=52+1n=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 ee. Soit en effet e me_m le plus petit exposant des flottants normalisés et considérons les flottants x=b e m(1+1/b)x=b^{e_m}(1+1/b) et y=b e my=b^{e_m}. Ces flottants sont distincts, mais leur différence n’est plus représentable par un flottant normalisé. Comme on ne souhaite pas représenter xyx-y par 0, (puisque le test x==yx==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,b n1[[0,b^{n-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/101/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/101/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=±1s=\pm 1 sur 1 bit, la mantisse M[0,2 52[M \in [0,2^{52}[ sur 52 bits, et l’exposant e[0,2 11[e \in [0, 2^{11}[ sur 11 bits. Pour les nombres “normaux”, l’exposant est en fait compris entre 1 et 2 1122^{11}-2, le nombre représenté est le rationnel (1+M2 52)2 e+12 10 (1+\frac{M}{2^{52}}) 2^{e+1-2^{10}} 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 rr dans [1,2[[1,2[, ce qui permet de déterminer l’exposant ee. Ensuite on écrit la représentation en base 2 de r1[0,1[r-1 \in [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=0e=0 ou e=2 111e=2^{11}-1) ont été introduites pour représenter ±\pm \infty (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

Voici un programme C++ affichant la représentation interne des flottants

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)f(x) avec xx mal connu de calculer aussi f(x)f'(x), en effet : f(x(1+h))=f(x)+xhf(x)+O(h 2)f(x(1+h))= f(x)+xh f'(x) + O(h^2) l’erreur relative sur f(x)f(x) est donc au premier ordre multipliée par |xf(x)f(x)||\frac{xf'(x)}{f(x)}| Par exemple, l’erreur relative sur e xe^x est au premier ordre l’erreur relative sur xx multipliée par |x||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 |xx 0|ε 0|x -x_0| \leq \varepsilon_0 et si |yy 0|ε 1|y-y_0| \leq \varepsilon_1 alors par l’inégalité triangulaire (|a+b||a|+|b||a+b|\leq |a|+|b|), on a : |(x+y)(x 0+y 0)||xx 0|+|yy 0|ε 0+ε 1 |(x+y)-(x_0+y_0)| \leq |x-x_0| + | y-y_0 | \leq \varepsilon_0 + \varepsilon_1 on dit que les erreurs absolues s’additionnent.

Définition 2   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 : |xx 0|ε |x-x_0| \leq \varepsilon

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

Définition 3   L’erreur relative est égale à l’erreur absolue divisée par la valeur absolue du nombre |xx 0|ε|x 0| |x-x_0| \leq \varepsilon |x_0|

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 532^{-53}, à titre d’exercice, on pourra vérifier que cette erreur d’arrondi est majorée par l’erreur absolue de la somme x+yx+y dès l’instant où xx et yy ont eux-même une erreur d’arrondi).

Lorsqu’on effectue une multiplication de deux nombres x,yx,y dont les représentants x 0,y 0x_0,y_0 sont non nuls, on a |xyx 0y 0x 0y 0|=|xx 0yy 01|=|(xx 01)(yy 01)+(xx 01)+(yy 01)| \left| \frac{xy-x_0 y_0}{x_0 y_0} \right| = \left| \frac{x}{x_0} \frac{y}{y_0} -1 \right| = \left| (\frac{x}{x_0}-1)(\frac{y}{y_0} -1)+(\frac{x}{x_0}-1)+(\frac{y}{y_0} -1) \right| 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 532^{-53} sur x 0y 0x_0 y_0.

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 xx et yy non nuls est de 2 532^{-53}, alors l’erreur relative sur xyxy sera de 2 53+2 53+2 106+2 533×2 53 2^{-53} + 2^{-53} + 2^{-106} + 2^{-53} \approx 3 \times 2^{-53} Lorsque l’erreur relative sur les données est grande devant 2 532^{-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 xx est représenté par x 0=1+2 52x_0=1+2^{-52} avec une erreur d’arrondi de 2 532^{-53} et yy par y 0=1y_0=-1 avec la même erreur d’arrondi, l’addition de xx et yy renvoie 2 522^{-52} avec une erreur absolue de 2*2 532 * 2^{-53} (ici il n’y a pas d’arrondi lorsqu’on fait la somme). C’est une erreur relative de 11 (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.01+2 52 (2.0^{-53}+2.0^{-53})+1.0 \rightarrow 1+2^{-52} alors que 2.0 53+(2.0 53+1.0)1 2.0^{-53}+(2.0^{-53}+1.0) \rightarrow 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)(f(x+h)-f(x-h))/(2h) que (f(x+h)f(x))/h(f(x+h)-f(x))/h car le terme d’erreur est en O(h 2)O(h^2) et non en O(h)O(h). Attention toutefois à ne pas prendre hh trop petit, sinon x+h=xx+h=x en flottants et même si x+hxx+h \neq x, l’erreur absolue sur f(x+h)f(xh)f(x+h)-f(x-h) est (au moins) d’ordre ε|f(x)|\varepsilon |f(x)|, donc l’erreur relative est d’ordre ε/h|f(x)|\varepsilon/h |f(x)|. Par exemple pour h=1e-8 le reste est en O(h 2)O(h^2) donc de l’ordre des erreurs d’arrondi mais l’erreur relative sur f(x+h)f(xh)f(x+h)-f(x-h) est d’ordre ε/h\epsilon/h largement supérieure (en flottants double-précision). On choisira plutôt hh tel que ε/h\epsilon/h soit proche de h 2h^2, 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))\exp(\sin(x)) en x=1x=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 nn 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 P(|X|>\alpha) \leq \frac{n\sigma^2}{\alpha^2} XX est la variable aléatoire somme des nn erreurs, α\alpha l’erreur et nσ 2n\sigma^2 la variance de la somme nn erreurs supposées indépendantes, cette probabilité tend vers 0 pour nn grand si α\alpha est d’ordre nn, et ne tend pas vers 0 si α\alpha est de l’ordre de n\sqrt{n}. Exemple : somme de n=400n=400 nombres répartis sur [1,1][-1,1] selon la loi uniforme (représentant des erreurs), on divise par n=\sqrt{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.

Attention, si on effectue la somme de nn réels jx j\sum_j x_j, les erreurs d’arrondis ne satisfont pas à ces hypothèses. En effet, l’erreur d’arrondi à chaque opération est une erreur relative, l’erreur absolue correspondante est ε|x 1+x 2|\epsilon |x_1+x_2| puis ε|x 1+x 2+x 3|\epsilon |x_1+x_2+x_3| puis ... ε|x 1+x 2+...+x n|\epsilon |x_1+x_2+...+x_n|, que l’on peut majorer par ε((n1)|x 1|+(n2)|x 2|+...+|x n||)\epsilon ((n-1)|x_1|+(n-2)|x_2|+...+|x_n||) La majoration de l’erreur d’arrondi dépend donc de l’ordre des termes, on a intérêt à sommer en commençant par les termes les plus petits en valeur absolue. Voici des programmes C++ illustrant cela

Mais on peut faire mieux, il est possible de corriger les erreurs d’arrondi dans une somme avec le programme suivant pour une liste (on peut bien sur adapter à la somme d’une expression dépendant d’une variable entière sans stocker de liste) :

Somme(l):={
  local x,s,c;
  s:=0.0;
  c:=0.0;
  pour x in l faire
    c += (x-((s+x)-s));
    s += x;
  fpour;
  print(c);
  return s+c;
}:;

onload
En effet, cc devrait valoir 0 sans erreurs d’arrondis, avec les erreurs d’arrondis, on a le premier calcul s+xs+x qui donnera une erreur opposée à celui du calcul de ss à la ligne suivante, le 2ième calcul effectué (s+x)s(s+x)-s donne une erreur absolue en ε|x|\epsilon |x| au pire (car c’est une erreur relative par rapport à (s+x)s(s+x)-s), et la 3ième erreur d’arrondi est négligeable (puisque la somme vaut 0). On a donc une erreur absolue sur s+cs+c qui est au premier ordre au pire en O(ε|x i|)O(\epsilon \sum|x_i|), bien meilleure que la majoration ε((n1)|x 1|+(n2)|x 2|+...+|x n||)\epsilon ((n-1)|x_1|+(n-2)|x_2|+...+|x_n||) calculée précédemment.

Par exemple
à comparer avec
(le calcul de SS est fait en exact, celui de sum(1. /j,j,1,n) est approché sans correction).

En conclusion, il est souvent très difficile de calculer une majoration rigoureuse de l’erreur pour des calculs (sauf les plus simples), et cette majoration est en général bien trop pessimiste. 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]xy, \ x \in [-1,2], y \in [-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 P=P_+-P_- en fonction du signe, puis utilise la monotonie de P +P_+ et P P_- sur +\mathbb{R}^+ et \mathbb{R}^- respectivement.

L’arithmétique d’intervalle dans \mathbb{C} 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]xy, \ x \in [-1,2], y \in [-1,2] et x 2,x[1,2]x^2, \ x \in [-1,2].

2.4  Types composés.

2.4.1  Les rationnels.

On sait donc représenter les entiers, pour les rationnels, il suffit de les représenter comme un couple d’entiers correspondant à leur écriture sous forme de fraction irréductible avec un dénominateur positif.

Proposition 4   L’algorithme d’Euclide permet de calculer le PGCD (plus grand commun diviseur) de 2 entiers, écrit ici en syntaxe Xcas :
pgcd(x,y):={
  local r;
  while (y!=0){
    r:=irem(x,y); // reste de x par y
    x:=y; // PGCD(x,y)=PGCD(y,r) donc on decale
    y:=r;
  }
  return x; // c'est le resultat car PGCD(x,0)=x
}

Preuve : on utilise le fait qu’un nombre divise aa et bb si et seulement si il divise r=abqr=a-bq et bb. Le PGCD de aa et bb est donc le PGCD de bb et du reste de la division euclidienne de aa par bb. Comme le reste est en valeur absolue plus petite que |b||b|, la taille des variables x,y,r décroit à chaque itération. Arrive un moment où le reste est nul, le PGCD est alors l’entier par lequel on a divisé. Il existe des variantes de cet algorithme un peu plus efficaces lorsque les nombres sont représentés en base 2 (PGCD binaire, voir par exemple A. Cohen).

On utilise cet algorithme et la division euclidienne pour simplifier une fraction d’entiers par le PGCD du numérateur et du dénominateur pour l’écrire sous forme irréductible.

Les calculs sont maintenant exacts et sans limitation de capacité (ou presque, la taille des entiers longs est bornée parce que la taille du champ mémoire fixant la longueur de stockage est bornée) mais souvent trop lents pour les calculs numériques usuels (par exemple pour calculer la valeur approchée de cosinus 23 degrés 27 minutes). On utilise alors un autre type dont les calculs de base sont souvent gérés par le microprocesseur (ou son coprocesseur arithmétique).

2.4.2  Les complexes

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.

2.4.3  Les polynômes

Après les nombres, l’objet le plus utilisé dans les systèmes de calcul formel est probablement le polynôme, toute simplification d’une expression se ramène à un moment donné à mettre sous forme irréductible une fraction de polynômes. Les principales représentations possibles sont :

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, il existe des algorithmes plus efficaces, cf. le chapitre sur les polynômes) Lien avec la représentation en base zz (TD).

Les polynômes peuvent servir à représenter des nombres non rationnels de manière exacte, par exemple les nombres algébriques (qui sont solutions d’une équation polynomiale).

2.4.4  Calcul symbolique

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.

2.4.5  Listes, séquences, tables

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

3.1  Rappel : suite récurrente

Une suite récurrente à un cran est définie par une valeur initiale u 0u_0 et une relation de récurrence que l’on va supposer ne pas dépendre explicitement de nn, donc de la forme u n+1=f(u n)u_{n+1}=f(u_n) On peut représenter graphiquement une suite de ce type par un graphe “en toile d’araignée”, par exemple pour les cinq premiers termes de f(x)=x+2x+1f(x)=\frac{x+2}{x+1} et u 0=3u_0=3,

gl_x=-0.5..2.5;plotseq((x+2)/(x+1),[0,0,2],5)

onload
La suite peut converger ou pas, par exemple avec f(x)=tan(x)0.1f(x)=\tan(x)-0.1 et f(x)=cos(x)f(x)=\cos(x)

3.2  Le point fixe

Soit ff une fonction continue sur un intervalle I=[a,b]I=[a,b] de \mathbb{R}, et à valeurs dans II (attention à bien choisir II pour que l’image de II par ff reste dans II). On s’intéresse à la suite u n+1=f(u n),u 0I(1) u_{n+1}=f(u_n), \quad u_0 \in I \qquad (1) Supposons que u nu_n converge vers une limite lIl \in I lorsque n+n \rightarrow +\infty, alors la limite doit vérifier f(l)=l f(l)=l puisque ff est continue. On dit que ll est un point fixe de ff. Ceci amène à l’idée d’utiliser ces suites pour résoudre numériquement l’équation f(x)=xf(x)=x. Nous allons donner un théorème permettant d’assurer que la suite (1) converge, et que la limite est l’unique solution de f(l)=lf(l)=l sur II.

Définition 1   On dit que ff est contractante de rapport k<1k&lt;1 sur II si ff est définie de II dans II et x,yI,|f(y)f(x)|k|yx| \forall x,y \in I, \quad |f(y)-f(x)| \leq k |y-x|

En pratique, les fonctions ff 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] f(y)-f(x)=f'(\theta) (y-x), \quad \theta \in [x,y] ainsi pour vérifier que ff est contractante, on étudie la valeur absolue de ff' sur II, il suffit de montrer que cette valeur absolue est strictement inférieure à un réel k<1k&lt;1 pour conclure (il faut donc chercher le maximum de |f||f'| sur II. Attention, il s’agit du maximum de |f||f'| et pas du maximum de ff', ce qui revient à chercher le maximum de ff' et de f-f').

On a alors le

Théorème 2   (du point fixe)
si
ff est contractante de I=[a,b]I=[a,b] dans II de rapport kk alors la suite (1) converge vers l’unique solution de f(l)=lf(l)=l dans II. On a de plus les encadrements : |u nl|k n|ba|,|u nl||u n+1u n|1k(2) |u_n-l| \leq k^n |b-a|, \quad |u_n -l | \leq \frac{|u_{n+1}-u_n|}{1-k} \qquad (2)

Démonstration : Tout d’abord si ff est contractante, on montre à partir de la définition de la continuité que ff est continue. Soit g(x)=f(x)xg(x)=f(x)-x, alors gg est continue, positive en aa et négative en bb, il existe donc l[a,b]l\in[a,b] tel que g(l)=0g(l)=0 (théorème des valeurs intermédiaires). Soit u nu_n une suite définie par (1). On a alors pour tout nn |u n+1l|=|f(u n)f(l)|k|u nl| |u_{n+1}-l|=|f(u_n)-f(l)| \leq k |u_n-l| Donc par une récurrence évidente : |u nl|k n|u 0l| |u_n-l| \leq k^n |u_0-l| ce qui entraine d’ailleurs que |u nl|k n|ab||u_n-l| \leq k^n |a-b|. Comme k[0,1[k \in [0,1[ , la suite géométrique k nk^n converge vers 0 lorsque nn tend vers l’infini, donc u nu_n tend vers ll. Notons que ll est unique car si ll' est une autre solution alors |ll|=|f(l)f(l)|k|ll||l-l'|=|f(l)-f(l')| \leq k|l-l'| donc (1k)|ll|0(1-k)|l-l'| \leq 0, or 1k>01-k&gt;0 et |ll|0|l-l'| \geq 0 donc |ll||l-l'| doit être nul.

Passons à la preuve de la majoration (2) 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>0m&gt;0 : u nl=u nu n+1+u n+1u n+2+...+u n+m1u n+m+u ml u_n-l= u_n - u_{n+1} + u_{n+1} - u_{n+2} + ... + u_{n+m-1}- u_{n+m} + u_m-l puis on majore avec l’inégalité triangulaire |u nl| j=0 m1|u n+ju n+j+1|+|u ml| |u_n-l| \leq \sum_{j=0}^{m-1} |u_{n+j}-u_{n+j+1}| + |u_m-l| puis on applique le fait que ff est contractante de rapport kk |u nl| j=0 m1k j|u nu n+1|+|u ml| |u_n-l| \leq \sum_{j=0}^{m-1} k^j |u_{n}-u_{n+1}| + |u_m-l| soit |u nl|1k m1k|u nu n+1|+|u ml| |u_n-l| \leq \frac{1-k^m}{1-k} |u_{n}-u_{n+1}| + |u_m-l| On fait alors tendre mm vers l’infini d’où le résultat.

Exemple 1 :
Cherchons une valeur approchée de 2\sqrt{2} par cette méthode. Il faut d’abord trouver une fonction ff dont 2\sqrt{2} est un point fixe, par exemple f(x)=x+2x+1 f(x)=\frac{x+2}{x+1} On vérifie que f(2)=2)f(\sqrt{2})=\sqrt{2}), puis que f=1/(x+1) 2f'=-1/(x+1)^2 donc ff décroit. On va voir si les hypothèses du théorème du point fixe s’appliquent sur par exemple [1,2][1,2]. Comme ff est décroissante f([1,2])=[f(2),f(1)]=[4/3,3/2]f([1,2])=[f(2),f(1)]=[4/3,3/2] qui est bien inclus dans [1,2][1,2] . De plus ff' est comprise entre 1/(1+1) 2=1/4-1/(1+1)^2=-1/4 et 1/(2+1) 2=1/9-1/(2+1)^2=-1/9 donc |f|<1/4|f'|&lt;1/4, ff est contractante de rapport 1/41/4. On peut donc itérer la suite à partir par exemple de u 0=1u_0=1 et on va converger vers 2\sqrt{2} (en s’en rapprochant à chaque cran d’un rapport inférieur à 1/41/4).



Exemple 2 :
Considérons l’équation en xx xesin(x)=t,e[0,1[ x- e \sin(x) =t, \quad e \in [0,1[ c’est l’équation du temps utilisée en astronomie pour trouver la position (plus précisément l’angle avec le grand axe de l’ellipse) d’une planète à l’instant tt sur son orbite elliptique (ee étant l’excentricité de l’ellipse).

gl_x=-3.5..3.5;gl_y:=-2.2..2.2; 
E:=ellipse(-1,1,2,color=red); cercle(0,2); theta:=2; 
line(0,slope=tan(theta)); angle(0,5,5*exp(i*theta),"x");
D:=droite(x=2*cos(theta)); M:=inter(D,E,i,legend="M(t)");
E:=ellipse(-1,1,2,color=red); cercle(0,2); theta:=2; 
line(0,slope=tan(theta)); angle(0,5,5*exp(i*theta),"x");
D:=droite(x=2*cos(theta)); M:=inter(D,E,i,legend="M(t)",color=red)

onload
Il n’y a pas de formule exacte permettant de calculer xx en fonction de tt. Si on a une valeur numérique pour tt, on peut trouver une valeur numérique approchée de xx par la méthode du point fixe, en réécrivant l’équation sous la forme f(x)=t+esin(x)=x f(x)=t+e\sin(x) = x On observe que ff envoie \mathbb{R} dans [te,t+e][t-e,t+e] donc on peut prendre I=[te,t+e]I=[t-e,t+e], de plus |f|e<1|f'|\leq e &lt;1, ff est contractante de rapport e[0,1[e \in [0,1[, le théorème s’applique, il suffit de prendre une valeur initiale dans [te,t+e][t-e,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 xx à 10 610^{-6} près, il suffira que la différence entre deux termes successifs de la suite u nu_n vérifie |u n+1u n|10 6(1e) |u_{n+1}-u_n| \leq 10^{-6} (1-e) on aura alors bien : |u nx||u n+1u n|1e10 6 |u_n-x| \leq \frac{|u_{n+1}-u_n|}{1-e} \leq 10^{-6} Par exemple, pour t=π/2t=\pi/2 et E=0.0167E=0.0167 (la Terre), au bout de 3 itérations

Cette méthode n’est pas toujours optimale, car la vitesse de convergence vers la limite ll est dite “linéaire”, c’est-à-dire que le temps de calcul pour avoir nn décimales est proportionnel à nn (ou encore il faut effectuer un nombre d’itérations proportionnel à nn, chaque itération faisant gagner en précision de l’ordre du rapport kk de contractance). En effet, supposons que ff' est continue en ll et que 0<L=|f(l)|<10&lt;L=|f'(l)|&lt;1 . Il existe alors un intervalle I=[lη,l+η]I=[l-\eta,l+\eta] tel que xIL2|f(x)|1+L2 x \in I \Rightarrow \frac{L}{2} \leq |f'(x)| \leq \frac{1+L}{2} Le théorème des accroissements finis donne alors |u n+1l|=|f(u n)f(l)|=|f(θ)||u nl|,θ[u n,l] |u_{n+1} - l | = |f(u_n)-f(l)| = |f'(\theta)| |u_n-l|, \quad \theta \in [u_n,l] Si u 0Iu_0 \in I, alors θI\theta \in I donc |u 1l||u 0l||u_1-l| \leq |u_0-l| et u 1Iu_1 \in I, par récurrence on a pour tout nn, u nIu_n \in I L2|u nl||u n+1l|1+L2|u nl| \frac{L}{2} |u_n-l| \leq |u_{n+1} - l| \leq \frac{1+L}{2} |u_n-l| on a donc par récurrence (L2) n|u 0l||u nl|(1+L2) n|u 0l| \left(\frac{L}{2}\right)^n|u_0-l| \leq |u_n-l| \leq \left( \frac{1+L}{2} \right)^n|u_0-l| Donc pour avoir |u nl|ε|u_n-l| \leq \epsilon il suffit que (1+L2) n|u 0l|εnln(ε|u 0l|)ln(1+L2) \left( \frac{1+L}{2} \right)^n|u_0-l| \leq \epsilon \Rightarrow n \geq \frac{\ln(\frac{\epsilon}{|u_0-l|})}{\ln( \frac{1+L}{2}) } et il faut que (L2) n|u 0l|εnln(ε|u 0l|)ln(L2) \left(\frac{L}{2}\right)^n |u_0-l| \leq \epsilon \Rightarrow n \geq \frac{\ln(\frac{\epsilon}{|u_0-l|})}{\ln( \frac{L}{2}) }

Si ff 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/21/2 (y compris prés de la solution de f(x)=xf(x)=x).

def dicho(f,a,b,eps):
    # local c,niter;
    if f(a)*f(b)>=0: 
        return "erreur: f(a)*f(b)>=0"
    while b-a>eps:
        c=(a+b)/2.0
        if f(a)*f(c)>0:
            a=c
        else:
            b=c
    return c

onload

Cette méthode a toutefois l’avantage de se généraliser en dimension supérieure à 1, contrairement à la dichotomie.

3.3  La méthode de Newton.

La méthode de Newton est une méthode de résolution de l’équation f(x)=0f(x)=0, attention à la différence avec le théorème du point fixe qui permet de résoudre numériquement f(x)=xf(x)=x. Si x 0x_0 est proche de la racine rr on peut faire un développement de Taylor à l’ordre 1 de la fonction ff en x 0x_0 : f(x)=f(x 0)+(xx 0)f(x 0)+O((xx 0) 2) f(x)=f(x_0)+(x-x_0)f'(x_0)+O((x-x_0)^2) Pour trouver une valeur approchée de rr, on ne garde que la partie linéaire du développement, on résout : f(r)=0f(x 0)+(rx 0)f(x 0) f(r)=0 \approx f(x_0) + (r-x_0) f'(x_0) donc (si f(x 0)0f'(x_0)\neq 0) : rx 0f(x 0)f(x 0) r \approx x_0 -\frac{f(x_0)}{f'(x_0)} Graphiquement, cela revient à tracer la tangente à la courbe représentative de ff et à chercher où elle coupe l’axe des xx.

gl_x=-3..3;
f:=x^2-2; 
G:=plot(f,x=-3..3);
x0:=1.9; M:=element(G,x0);
T:=tangent(G,x0,legend="",color=red);
N:=inter_unique(T,droite(y=0));
legend(-3,"xN="+round(abscisse(N),3))

onload
On considère donc la suite récurrente définie par une valeur u 0u_0 proche de la racine et par la relation : u n+1=u nf(u n)f(u n) u_{n+1} = u_n -\frac{f(u_n)}{f'(u_n)} Remarque : C’est un cas particulier de point fixe pour l’équation g(x)=xg(x)=x avec g(x)=xf(x)f(x),g(x)=f(x)f(x)f(x) 2g(x)=x-\frac{f(x)}{f'(x)}, \quad g'(x)=-\frac{f(x) f'{'}(x)}{f'(x)^2} et donc gg' s’annule en un point où ff est nulle. On espère ainsi avoir une convergence très rapide près d’une racine.

Il y a deux théorèmes importants, l’un d’eux prouve que si u 0u_0 est “assez proche” de rr alors la suite u nu_n converge vers rr, malheureusement il est difficile de savoir en pratique si on est “assez proche” de u 0u_0 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 3   Soit ff une fonction de classe C 2C^2 (2 fois continument dérivable) sur un intervalle fermé II. Soit rr une racine simple de ff située à l’intérieur de II (telle que f(r)=0f(r)=0 et f(r)0f'(r)\neq 0). Alors il existe ε>0\varepsilon&gt;0 tel que la suite définie par u n+1=u nf(u n)f(u n),|u 0r|ε u_{n+1} = u_n -\frac{f(u_n)}{f'(u_n)}, \quad |u_0-r| \leq \varepsilon converge vers rr.

Si on a |f|M|f'{'}| \leq M et |1/f|m|1/f'| \leq m sur un intervalle [rη,r+η][r-\eta,r+\eta] contenu dans II, alors on peut prendre tout réel ε>0\varepsilon&gt;0 tel que ε<2/(mM)\varepsilon &lt; 2/(mM) et εη\varepsilon \leq \eta.

Démonstration : on a u n+1r=u nrf(u n)f(u n)=(u nr)f(u n)f(u n)f(u n) u_{n+1}-r = u_n - r - \frac{f(u_n)}{f'(u_n)} = \frac{(u_n-r)f'(u_n)-f(u_n)}{f'(u_n)} En appliquant un développement de Taylor de ff en u nu_n à l’ordre 2, on obtient pour un réel θ\theta situé entre rr et u nu_n : 0=f(r)=f(u n)+(ru n)f(u n)+(ru n) 2f(θ)2 0 = f(r)=f(u_n)+(r-u_n) f'(u_n) + (r-u_n)^2 \frac{f'{'}(\theta)}{2} donc : (u nr)f(u n)f(u n)=(u nr) 2f(θ)2 (u_n-r)f'(u_n)-f(u_n)= (u_n-r)^2 \frac{f'{'}(\theta)}{2} d’où : |u n+1r||u nr| 21|f(u n)||f(θ)|2 |u_{n+1}-r| \leq |u_n-r|^2 \frac{1}{|f'(u_n)|} \frac{|f'{'}(\theta)|}{2} On commence par choisir un intervalle [rε,r+ε][r-\varepsilon,r+\varepsilon] contenant strictement rr et tel que |f|<M|f'{'}|&lt;M et |1/f|<m|1/f'|&lt;m sur [rε,r+ε][r-\varepsilon,r+\varepsilon] (c’est toujours possible car ff'{'} et 1/f1/f' sont continues au voisinage de rr puisque f(r)0f'(r)\neq 0). Si u nu_n est dans cet intervalle, alors θ\theta aussi donc |u n+1r||u nr| 2Mm2|u nr|Mm2|u nr|, |u_{n+1}-r| \leq |u_n-r|^2 \frac{Mm}{2} \leq \frac{|u_n-r|Mm}{2} |u_n-r|, On a |u nr|ε|u_n-r| \leq \varepsilon, on diminue si nécessaire ε\varepsilon pour avoir ε<2/(Mm)\varepsilon &lt; 2/(Mm), on a alors : |u n+1r|k|u nr|,k=εMm2<1 |u_{n+1}-r| \leq k |u_n-r|, \quad k=\frac{\varepsilon Mm}{2}&lt;1 donc d’une part u n+1u_{n+1} est encore dans l’intervalle [rε,r+ε][r-\varepsilon,r+\varepsilon] 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 rr. En fait la convergence est bien meilleure lorsqu’on est proche de rr grace au carré dans |u nr| 2|u_n-r|^2, plus précisément, on montre par récurrence que |u nr||u 0r| 2 n(Mm2) 2 n1 |u_n-r| \leq |u_0 - r|^{2^n} \left( \frac{Mm}{2} \right)^{2^n-1} il faut donc un nombre d’itérations proportionnel à ln(n)\ln(n) pour atteindre une précision donnée.

Remarque : ce théorème se généralise sur \mathbb{C} et même sur n\mathbb{R}^n.

Exemple : pour calculer 2\sqrt{2}, on écrit l’équation x 22=0x^2-2=0 qui a 2\sqrt{2} comme racine simple sur I=[1/2,2]I=[1/2,2], on obtient la suite récurrente u n+1=u nu n 222u n u_{n+1} = u_n - \frac{u_n^2-2}{2u_n} Si on prend η=1/2\eta=1/2, on a f=2xf'=2x et f=2f'{'}=2 donc on peut prendre M=2M=2 et m=1m=1 car |1/f|1|1/f'|\leq 1 sur [21/2,2+1/2][\sqrt{2}-1/2,\sqrt{2}+1/2]. On a 2/(mM)=12/(mM)=1, on peut donc prendre ε=1/2\varepsilon=1/2, la suite convergera pour tout u 0[21/2,2+1/2] u_0 \in [\sqrt{2}-1/2,\sqrt{2}+1/2].

Plus générallement, on peut calculer une racine kk-ième d’un réel aa en résolvant f(x)=x kaf(x)=x^k-a 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+1n+1 couleurs selon que la suite définie par la méthode de Newton converge vers l’une des nn racines d’un polynôme de degré nn fixé au bout de par exemple 50 itérations (la n+1n+1-ième couleur servant aux origines de suite qui ne semblent pas converger). Session Xcas

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

Définition 4   (convexité)
Une fonction
ff continument dérivable sur un intervalle II de \mathbb{R} est dite convexe si son graphe est au-dessus de la tangente en tout point de II.

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

Théorème 5   Si ff est C 2C^2 et f0f'{'} \geq 0 sur II alors ff est convexe.

Démonstration :
L’équation de la tangente au graphe en x 0x_0 est y=f(x 0)+f(x 0)(xx 0) y=f(x_0)+f'(x_0)(x-x_0) Soit g(x)=f(x)(f(x 0)+f(x 0)(xx 0)) g(x)=f(x)-(f(x_0)+f'(x_0)(x-x_0)) on a : g(x 0)=0,g(x)=f(x)f(x 0),g(x 0)=0,g=f0 g(x_0)=0, \quad g'(x)=f'(x)-f'(x_0), \quad g'(x_0)=0, \quad g'{'}=f'{'} \geq 0 donc gg' est croissante, comme g(x 0)=0g'(x_0)=0, gg' est négative pour x<x 0x&lt;x_0 et positive pour x>x 0x&gt;x_0, donc gg est décroissante pour x<x 0x&lt;x_0 et croissante pour x>x 0x&gt;x_0. On conclut alors que g0g \geq 0 puisque g(x 0)=0g(x_0)=0. Donc ff est bien au-dessus de sa tangente.

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

Théorème 6   Si f(r)=0,f(r)>0f(r)=0, f'(r)&gt;0 et si f0f'{'} \geq 0 sur [r,b][r,b] alors pour tout u 0[r,b]u_0 \in [r,b] la suite de la méthode de Newton u n+1=u nf(u n)f(u n), u_{n+1} = u_n -\frac{f(u_n)}{f'(u_n)}, est définie, décroissante, minorée par rr et converge vers rr. De plus 0u nrf(u n)f(r) 0 \leq u_n -r \leq \frac{f(u_n)}{f'(r)}

Démonstration :
On a f0f'{'} \geq 0 donc si f(r)>0f'(r)&gt;0 alors f>0f'&gt;0 sur [r,b][r,b], ff est donc strictement croissante sur [r,b][r,b] on en déduit que f>0f&gt;0 sur ]r,b]]r,b] donc u n+1u nu_{n+1} \leq u_n. Comme la courbe représentative de ff est au-dessus de la tangente, on a u n+1ru_{n+1} \geq r (car u n+1u_{n+1} est l’abscisse du point d’intersection de la tangente avec l’axe des xx). La suite u nu_n est donc décroissante minorée par rr, donc convergente vers une limite lrl\geq r. À la limite, on a l=lf(l)f(l)f(l)=0 l=l-\frac{f(l)}{f'(l)} \Rightarrow f(l)=0 donc l=rl=r car f>0f&gt;0 sur ]r,b]]r,b].

Comme (u n)(u_n) est décroissante, on a bien 0u nr0 \leq u_n -r, pour montrer l’autre inégalité, on applique le théorème des accroissements finis, il existe θ[r,u n]\theta \in [r,u_n] tel que f(u n)f(r)=(u nr)f(θ) f(u_n)-f(r)=(u_n-r)f'(\theta) comme f(r)=0f(r)=0, on a u nr=f(u n)f(θ) u_n-r = \frac{f(u_n)}{f'(\theta)} et la deuxième inégalité du théorème en découle parce que ff' est croissante.

Variantes :
Il existe des variantes, par exemple si f(r)<0f'(r)&lt;0 et f0f'{'} \geq 0 sur [a,r][a,r]. Si f0f'{'} \leq 0, on considère g=fg=-f.

Application :
On peut calculer la valeur approchée de la racine kk-ième d’un réel a>0a&gt;0 en appliquant ce deuxième théorème. En effet si a>0a&gt;0, alors x kax^k-a est 2 fois continument dérivable et de dérivée première kx k1kx^{k-1} et seconde k(k1)x k2k(k-1)x^{k-2} strictement positives sur R +*\R^{+\ast} (car k2k \geq 2). Il suffit donc de prendre une valeur de départ u 0u_0 plus grande que la racine kk-ième, par exemple 1+a/k1+a/k (en effet (1+a/k) k1+ka/k=1+a(1+a/k)^k \geq 1+k a/k=1+a). En appliquant l’inégalité du théorème, on a : 0u naku n kakak k1u n kakaaku n kaka(1+ak) 0 \leq u_n - \sqrt[k]{a} \leq \frac{u_n^k - a}{k\sqrt[k]{a}^{k-1} } \leq \frac{u_n^k-a}{ka} \sqrt[k]{a} \leq \frac{u_n^k-a}{ka} (1+\frac{a}{k}) Pour avoir une valeur approchée de ak\sqrt[k]{a} à ε\varepsilon près, on peut donc choisir comme test d’arrêt u n kaka1+akε u_n^k -a \leq \frac{ka}{1+\frac{a}{k}} \varepsilon Par exemple pour 2\sqrt{2}, le test d’arrêt serait u n 222εu_n^2-2 \leq 2 \varepsilon.

Remarque : Lorsque u nu_n est proche d’une racine, le calcul approché de f(u n)f(u_n) peut poser des problèmes à cause de compensations entre deux valeurs proches que l’on soustrait. Par exemple pour f(x)=x 22f(x)=x^2-2, si xx est proche de 2\sqrt{2}, x 2x^2 sera proche de 2, et le calcul de la différence va perdre énormément de précision relative.

On peut contourner ce problème en utilisant des flottants multipŕecision, ci-dessus en posant par exemple u=evalf(u,40). Dans le cas de la recherche de racines de polynômes, on peut utiliser des rationnels (u=2) mais un autre problème se produit : la taille du numérateur et du dénominateur va en gros doubler à chaque itération! Pour éviter cela, on peut remplacer par un rationnel proche de taille de numérateur/dénominateur plus petite par arrondi à une puissance de 2 dépendant de la précision attendue de u nu_n.

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

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

Soit ff une fonction indéfiniment dérivable sur un intervalle II de \mathbb{R} et x 0Ix_0 \in I. On peut alors effectuer le développement de Taylor de ff en x 0x_0 à l’ordre nn T n(f)(x)=f(x 0)+(xx 0)f(x 0)+...+(xx 0) nf [n](x 0)n! T_n(f)(x)= f(x_0) + (x-x_0) f'(x_0) + ... + (x-x_0)^n \frac{f^{[n]}(x_0)}{n!} et se demander si T n(f)T_n(f) converge lorsque nn tend vers l’infini, si la limite est égale à f(x)f(x) et si on peut facilement majorer la différence entre f(x)f(x) et T n(f)(x)T_n(f)(x). Si c’est le cas, on pourra utiliser T n(f)(x)T_n(f)(x) comme valeur approchée de f(x)f(x).

On peut parfois répondre à ces questions simultanément en regardant le développement de Taylor de ff avec reste 

Théorème 1   Si ff est de classe C n+1C^{n+1} sur [x 0,x][x_0,x] (ou [x,x 0][x,x_0]), il existe θ\theta compris entre x 0x_0 et xx tel que R n(x):=f(x)T n(f)(x)=(xx 0) n+1f [n+1](θ)(n+1)! R_n(x) := f(x)- T_n(f)(x) = (x-x_0)^{n+1}\frac{f^{[n+1]}(\theta)}{(n+1)!}

La preuve se fait en considérant la fonction g(t):=f(x)(f(t)+(xt)f(t)+...+(xt) nf [n](t)n!)C(xt) n+1g(t) := f(x)- (f(t) + (x-t) f'(t) + ... + (x-t)^n \frac{f^{[n]}(t)}{n!} ) - C(x-t)^{n+1} qui vérifie g(x)=0g(x)=0. On choisit CC pour que g(x 0)=0g(x_0)=0 et on calcule la dérivée qui s’annule en θ\theta entre x 0x_0 et xx, ce qui donne la valeur de CC.

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

4.1  La fonction exponentielle

Soit f(x)=exp(x)f(x)=\exp(x) et x 0=0x_0=0, la dérivée nn-ième de ff est exp(x)\exp(x), donc R n(x)=exp(θ)x n+1/(n+1)!R_n(x)=\exp(\theta)x^{n+1}/(n+1)! avec θ\theta compris entre 0 et xx, ainsi si xx est positif |R n(x)|e xx n+1/(n+1)!|R_n(x)| \leq e^x x^{n+1}/(n+1)! et si xx est négatif, |R n(x)|x n+1/(n+1)!|R_n(x)| \leq x^{n+1}/(n+1)!. Dans les deux cas, la limite de R nR_n est 0 lorsque nn tend vers l’infini, car pour n2xn \geq 2x, on a x n+1(n+1)!=x nn!xn+112x nn! \frac{x^{n+1}}{(n+1)!} = \frac{x^n}{n!} \frac{x}{n+1}\leq \frac{1}{2}\frac{x^n}{n!} on a donc pour tout xx réel e x=lim n+T n(f)(x)=lim n+ k=0 nx kk!= k=0 x kk! e^x = \lim_{n \rightarrow +\infty} T_n(f)(x) = \lim_{n \rightarrow +\infty} \sum_{k=0}^n \frac{x^k}{k!} = \sum_{k=0}^\infty \frac{x^k}{k!}

gl_x=-1.5..1.5,gl_y=-4..4
,plot([exp(x),1+x,1+x+x^2/2,1+x+x^2/2+x^3/6,
1+x+x^2/2+x^3/6+x^4/24],x=-1.5..1.5,
color=[black,blue,cyan,green,red])

onload

Comment en déduire une valeur approchée de e xe^x? Il suffira d’arrêter la sommation lorsque R:=(x) n+1/(n+1)!R:=(-x)^{n+1}/(n+1)! si x<0x&lt;0 ou lorsque R:=e xx n+1/(n+1)!R:=e^x x^{n+1}/(n+1)! si x>0x&gt;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 e xe^x à une erreur relative ε\varepsilon donnée (par exemple ε=2 53\varepsilon=2^{-53} pour stocker le résultat dans un double) il suffit que R/e x<εR/e^x &lt; \varepsilon, donc si xx est positif, il suffit que x n+1/(n+1)!<εx^{n+1}/(n+1)!&lt;\varepsilon, on peut donc arrêter la sommation lorsque le terme suivant est plus petit que ε\varepsilon.

On observe que plus xx est grand, plus nn 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 10e^{-10} ou plus générallement 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 532^{-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 e x=1/e xe^{-x}=1/e^x (ne change pas l’erreur relative). Pour les grands réels, on peut utiliser e 2x=(e x) 2e^{2x}=(e^x)^2 (multiplie par 2 l’erreur relative). On peut aussi, si on connait une valeur approchée de ln(2)\ln(2), effectuer la division euclidienne de xx par ln(2)\ln(2) avec reste symétrique : x=aln(2)+r,a,|r|ln(2)2 x = a \ln(2) + r, \quad a \in \mathbb{Z}, |r| \leq \frac{\ln(2)}{2} puis si rr est positif, on somme la série de T(f)(r)T(f)(r), si rr est négatif, on calcule T(f)(r)T(f)(-r) et on inverse, on applique alors : e x=2 ae r e^x = 2^a e^r

Il faut toutefois noter que ln(2)\ln(2) n’étant pas connu exactement, on commet une erreur d’arrondi absolu sur rr d’ordre aηa \eta, où η\eta est l’erreur relative sur ln(2)\ln(2), il faut donc ajouter une erreur d’arrondi relative de xη/ln(2)x \eta /\ln(2) qui peut devenir grande si xx est grand. Puis il faut ajouter la somme des erreurs d’arrondi due au calcul de e re^r, que l’on peut minimiser en utilisant la méthode de Horner pour évaluer T n(f)(r)T_n(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)\ln(2) avec cette précision, le calcul de e xe^x 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 532^{-53}).

Notons que en général xx lui-même a déjà été arrondi ou n’est connu qu’avec une précision relative. Or si x>0x&gt;0 est connu avec une erreur relative de ε\varepsilon (donc une erreur absolue de ε|x|\varepsilon |x|, alors e x+ε|x|=e xe ε|x| e^{x+\varepsilon |x|}= e^x e^{\varepsilon |x|} donc on ne peut pas espérer mieux qu’une erreur relative de e ε|x|1e^{\varepsilon |x|}-1 sur l’exponentielle de xx. Si εx\varepsilon x est petit cette erreur relative (impossible à éviter, quel que soit l’algorithme utilisé pour calculer l’exponentielle) est d’ordre ε|x|\varepsilon |x|. Si εx\varepsilon 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 xx 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+\pi)=-sin(x), sin(x)=sin(x)\sin(-x)=-\sin(x), sin(x)=cos(π/2x)\sin(x)=\cos(\pi/2-x) pour se ramener au calcul de sin(x)\sin(x) ou de cos(x)\cos(x) sur [0,π/4][0,\pi/4]. sin(x)= n=0 (1) nx 2n+1(2n+1)!,cos(x)= n=0 (1) nx 2n(2n)! \sin(x)=\sum_{n=0}^\infty (-1)^n \frac{x^{2n+1}}{(2n+1)!}, \quad \cos(x)=\sum_{n=0}^\infty (-1)^n \frac{x^{2n}}{(2n)!}

Cette méthode a toutefois ces limites, car il peut devenir impraticable de calculer la dérivée nn-ième d’une fonction (par exemple avec tan(x)\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.

Remarque : pour calculer la fonction exponentielle en optimisant le nombre d’opérations à effectuer, on utilise des approximants de Padé, ce sont des fractions rationnelles qui ont le développement de Taylor souhaité en x=0x=0, et présentent des symétries permettant de les calculer efficacement.

Le programme suivant permet de calculer l’exponentielle sur [ln(2)/2,ln(2)/2][-\ln(2)/2,\ln(2)/2] en 12 opérations avec une précision proche de 1e-15

def f(x): 
    X=x*x
    q=((X+420)*X+15120)*x
    r=(30X+3360)*X+30240
    return (r+q)/(r-q)

onload

En ajoutant une division par ln(2)\ln(2) avec arrondi à l’entier le plus proche, on a l’exponentielle en 16 opérations arithmétiques de base.

4.2  Séries entières.

Les séries de type prendre la limite lorsque nn tend vers l’infini du développement de Taylor en x=0 sont de la forme n=0 a nx n:=lim k+ n=0 ka nx n,a n=f [n](0)n! \sum_{n=0}^\infty a_n x^n := \lim_{ k \rightarrow +\infty} \sum_{n=0}^k a_n x^n, a_n=\frac{f^{[n]}(0)}{n!} On peut s’intéresser plus générallement à n=0 a nx n\sum_{n=0}^\infty a_n x^n lorsque a na_n 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 x 0x_0 tel que |a nx 0 n||a_n x_0^n| est borné (ce sera le cas en particulier si la série converge en x 0x_0), alors |a nx n|=|a nx 0 n||xx 0| nM|xx 0| n |a_n x^n| = |a_n x_0^n| |\frac{x}{x_0}|^n \leq M |\frac{x}{x_0}|^n la série converge donc en xx si |x|<|x 0||x|&lt;|x_0| et on peut majorer le reste de la série au rang nn par |R n|M|xx 0| n+11|xx 0| |R_n| \leq M \frac{ |\frac{x}{x_0}|^{n+1}} {1-|\frac{x}{x_0}|} 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 kk décimales dépend linéairement kk, les constantes sont d’autant plus grandes que |x||x| est grand).

Théorème 2   S’il existe un rang n 0n_0, un réel M>0M&gt;0 et un complexe x 0x_0 tels que pour n>n 0n&gt;n_0, on ait : |a nx 0| nM |a_n x_0|^n \leq M alors la série converge pour |x|<|x 0||x|&lt;|x_0| et pour nn 0n\geq n_0, on a : |R n|M|xx 0| n+11|xx 0|(3) |R_n| \leq M \frac{ |\frac{x}{x_0}|^{n+1}} {1-|\frac{x}{x_0}|} \qquad (3)

On en déduit qu’il existe un réel positif R0R\geq 0 éventuellement égal à ++\infty tel que la série converge (la limite de la somme jusqu’à l’infini existe) lorsque |x|<R|x|&lt;R et n’existe pas lorsque |x|>R|x|&gt;R, ce réel est appelé rayon de convergence de la série. Par exemple ce rayon vaut ++\infty pour l’exponentielle, le sinus ou le cosinus. Il est égal à 1 pour la série géométrique x n\sum x^n (car elle diverge si |x|>1|x|&gt;1 et converge si |x|<1|x|&lt;1). On ne peut pas dire ce qui se passe génériquement lorsqu’on est à la limite, c’est-à-dire lorsque |x|=R|x|=R (si R+R\neq +\infty). 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+R\neq +\infty est en gros la même que celle d’une série géométrique de raison |x|/R|x|/R.

Lorsque deux 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=1x+x 2x 3+... (1+x)^{-1} = 1-x+x^2-x^3+... et on obtient une série entière de rayon de convergence non nul. On peut aussi composer deux séries entières gg et ff en gfg\circ f (avec les règles de calcul de composition des polynômes) si f(0)=0f(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 a na_n est assez complexe. La fonction (1+x) α(1+x)^\alpha est développable en séries entières pour tout α\alpha \in \mathbb{R} avec un rayon de convergence 1 (ou l’infini pour α\alpha entier positif). (1+x) α=1+αx+α(α1)2!x 2+...+α(α1)...(αn+1)n!x n+... (1+x)^\alpha = 1 + \alpha x + \frac{\alpha (\alpha-1)}{2!} x^2 + ... + \frac{\alpha (\alpha-1) ... (\alpha -n +1)}{n!} x^n + ... Pour α=1\alpha=-1, c’est la série géométrique de raison x-x, en effet si |x|<1|x|&lt;1 : n=0 k(x) n=1(x) k+11+x k11+x \sum_{n=0}^k (-x)^n = \frac{1-(-x)^{k+1}}{1+x} \rightarrow_{k\rightarrow \infty} \frac{1}{1+x} En intégrant par rapport à xx, on obtient que ln(1+x)\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+1n+1 \ln(1+x) = \sum_{n=0}^\infty \frac{(-x)^{n+1}}{n+1} On peut calculer de manière analogue le développement en série entière de arctan(x)\arctan(x) en iintégrant celui de 1/(1+x 2)1/(1+x^2), de même pour arccos(x)\arccos(x) et arcsin(x)\arcsin(x) en intégrant celui de (1x 2) 1/2(1-x^2)^{-1/2}. arctan(x)= n=0 (1) nx 2n+12n+1, \arctan(x)=\sum_{n=0}^\infty (-1)^{n} \frac{x^{2n+1}}{2n+1}, On peut donc calculer ln\ln, arctan\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)\ln(1+x), on pourrait répondre comme avec l’exponentielle en majorant la dérivée n+1n+1-ième, mais ce n’est plus faisable pour arctan,arcsin,arccos\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.

4.3  Série alternée

Théorème 3   Soit S n= k=0 n(1) ku kS_n= \sum_{k=0}^n (-1)^k u_k la somme jusqu’au rang nn d’une série de réels tels que la suite des u ku_k décroit à partir d’un rang n 0n_0 et tend vers 0 lorsque k+k\rightarrow +\infty. Alors S nS_n converge vers une limite SS. Si nn 0n\geq n_0, la limite est comprise entre deux sommes partielles succesives S nS_n et S n+1S_{n+1} et le reste est majoré par la valeur absolue du premier terme non sommé : |R n||u n+1| |R_n| \leq |u_{n+1}|

Démonstration :
on montre que les suites v n=S 2nv_n=S_{2n} et w n=S 2n+1w_n=S_{2n+1} sont adjacentes. On a v n+1v n=S 2n+2S 2n=(1) 2n+2u 2n+2+(1) 2n+1u 2n+1=u 2n+2u 2n+10 v_{n+1}-v_n= S_{2n+2}-S_{2n}= (-1)^{2n+2} u_{2n+2} + (-1)^{2n+1} u_{2n+1} = u_{2n+2}-u_{2n+1} \leq 0 donc v nv_n est décroissante, de même w nw_n est croissante, et v nw n=u 2n+1v_n-w_n=u_{2n+1} est positif et tend vers 0. On en déduit que v nv_n et w nw_n convergent vers la même limite SS telle que v n>S>w nv_n&gt;S&gt;w_n 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 u nu_n tende assez vite vers 0, sinon il y aura perte de précision sur la mantisse lorsqu’on effectuera u 2nu 2n+1u_{2n}-u_{2n+1}. On sommera aussi les termes par ordre décroissant pour diminuer les erreurs d’arrondi.

4.4  La fonction logarithme

Si nous voulons calculer ln(1+x)\ln(1+x) pour x[0,1[x \in [0,1[ avec une précision ε\varepsilon, il suffit de calculer k=0 n(1) kx k+1k+1 \sum_{k=0}^n (-1)^k \frac{x^{k+1}}{k+1} pour nn tel que la valeur absolue du terme suivant soit plus petit que ε\varepsilon : ntel quex n+1n+1<ε n \ \mbox{tel que} \ \frac{x^{n+1}}{n+1} &lt; \varepsilon en effet, les signes sont alternés et la suite x k+1k+1\frac{x^{k+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 xx est voisin de 1, il faut calculer nn termes pour avoir une précision en 1/n1/n, par exemple 1 million de termes pour avoir une précision de 1e61e-6 (sans tenir compte des erreurs d’arrondi). Si xx est proche de 1/21/2 il faut de l’ordre de ln(ε)/ln(2)-\ln(\varepsilon)/\ln(2) termes ce qui est mieux, mais encore relativement grand (par exemple 50 termes environ pour une précision en 1e161e-16, 13 termes pour 1e41e-4). On a donc intérêt à se ramener si possible à calculer la fonction en un xx où la convergence est plus rapide (donc |x||x| le plus petit possible). Par exemple pour le calcul de ln(1+x)\ln(1+x) on peut :

Nous sommes donc en mesure de calculer précisément le logarithme ln(1+x)\ln(1+x) pour |x|<1/2|x|&lt;1/2. Pour calculer lnln sur +\mathbb{R}^+, on se ramène à [1/2,1][1/2,1] en utilisant l’écriture mantisse-exposant.

Remarquons que si xx est connu à une erreur relative ε\varepsilon près, comme ln(x(1±ε))=ln(x)+ln(1±ε) \ln(x(1 \pm \varepsilon))=\ln(x) + \ln(1 \pm \varepsilon) ln(x)\ln(x) est connu à une erreur absolue de |ln(1±ε)|ε|\ln(1 \pm \varepsilon)| \approx \varepsilon. Si ln(x)\ln(x) est proche de 0, on a une grande perte de précision relative.

Finalement, nous savons calculer ln\ln et exp\exp sous réserve d’avoir dans une table la valeur de ln(2)\ln(2). Pour calculer ln(2)\ln(2) précisément, on peut utiliser ln(2)=ln(1/2)=ln(11/2) \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)\ln(1+x) est borné par M=1M=1 en x=1x=1, donc d’après (3) : |R n|12 n |R_n| \leq \frac{1}{2^n} (on peut même obtenir 1/(n2 n)1/(n2^n) car on a besoin de MM uniquement pour les termes d’ordre plus grand que nn, on peut donc prendre M=1/nM=1/n). Par exemple, pour avoir ln(2)\ln(2) avec une mantisse de 80 bits, on effectue une fois pour toutes avec un logiciel de calcul formel :


puis la division en base 2 avec 81 bits de précision

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

4.5  Autres applications

On peut calculer certaines intégrales de la même manière, par exemple 0 1/211+x 3 \int _0^{1/2} \frac{1}{\sqrt{1+x^3}} mais aussi des fonctions définies par des intégrales (cas de nombreuses fonctions spéciales).

4.5.1  Exemple : la fonction d’erreur (error fonction, erf)

Cette fonction est définie à une constante multiplicative près par : f(x)= 0 xe t 2dt f(x)=\int_0^x e^{-t^2} \ dt On peut développer en séries entières l’intégrand (rayon de convergence ++\infty), puis intégrer terme à terme, on obtient f(x)= n=0 +(1) nx 2n+1n!(2n+1) f(x)= \sum_{n=0}^{+\infty} (-1)^n \frac{x^{2n+1}}{n! (2n+1)} Ce développement converge très rapidement pour |x|1|x|\leq 1. Par contre, pour |x||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)f(1) à 1e161e-16 près. Combien de termes faut-il calculer dans la somme pour trouver une valeur approchée de f(7)f(7) à 1e161e-16 près ? Comparer la valeur de f(7)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)f(7) par ce développement. Combien de termes faut-il calculer pour déterminer f(10)f(10) à 1e161e-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)f(10) ?

4.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 a n+1a_{n+1} et a na_n. Si on arrive à montrer par exemple qu’il y a une solution ayant un développement alternée, ou plus générallement, si on a une majoration |a n+1/a n|<C|a_{n+1}/a_n|&lt;C, alors le reste de la série entière est majoré par |a nx n|/(1|Cx|)|a_nx^n|/(1-|Cx|) lorsque |x|<1/C|x|&lt;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.

4.5.3  Exemple : fonctions de Bessel d’ordre entier

Soit mm un entier positif fixé, on considère l’équation différentielle x 2y+xy+(x 2m 2)y=0 x^2 y'{'} + x y' + (x^2-m^2)y=0 dont on cherche une solution série entière y= k=0 a kx ky=\sum_{k=0}^\infty a_k x^k . En remplacant dans l’équation, si xx est dans le rayon de convergence de la série (rayon supposé non nul), on obtient k=0 k(k1)a kx k+ k=0 ka kx k+ k=0 (x 2m 2)a kx k=0 \sum_{k=0}^\infty k(k-1)a_k x^k + \sum_{k=0}^\infty k a_k x^k + \sum_{k=0}^\infty (x^2-m^2) a_k x^k =0 soit encore 0 = k=0 (k 2m 2+x 2)a kx k = m 2a 0+(1m 2)a 1x+ k=2 [(k 2m 2)a k+a k2]x k \begin{matrix} 0 &=& \sum_{k=0}^\infty (k^2-m^2+x^2) a_k x^k \\ &=& -m^2 a_0 + (1-m^2)a_1 x + \sum_{k=2}^\infty [(k^2-m^2) a_k +a_{k-2}]x^k \end{matrix} Par exemple, prenons le cas m=0m=0. On a alors a 0a_0 quelconque, a 1a_1 nul et pour k2k\geq 2 a k=a k2k 2 a_k = - \frac{a_{k-2}}{k^2} Donc tous les aa d’indice impair sont nuls. Les pairs sont non nuls si a 00a_0\neq 0, et ils sont de signe alterné. Soit xx fixé, on observe que pour 2k>|x|2k &gt; |x|, |a 2kx 2k|<|a 2k2x 2k2| |a_{2k} x^{2k}| &lt; |a_{2k-2} x^{2k-2}| donc la série k=0 a kx k\sum_{k=0}^\infty a_k x^k est alternée à partir du rang partie entière de |x||x| plus un. Donc elle converge pour tout xx (le rayon de convergence de yy est ++\infty) et le reste de la somme jusqu’à l’ordre 2n2n est inférieur en valeur absolue à : |R 2n(x)||a 2n+2x 2n+2| |R_{2n}(x)| \leq |a_{2n+2} x^{2n+2}| Par exemple, pour avoir une valeur approchée à 1e101e-10 près de y(x)y(x) pour a 0=1a_0=1 et |x|1|x|\leq 1, on calcule y= k=0 2na kx ky=\sum_{k=0}^{2n} a_k x^k , on s’arrête au rang nn tel que |a 2n+2x 2n+2||a 2n+2|10 10 |a_{2n+2} x^{2n+2}| \leq |a_{2n+2}| \leq 10^{-10} On remarque que : a 2n=(1) n2 24 2...(2n) 2=(1) n2 2nn! 2 a_{2n} = \frac{(-1)^n}{2^2 4^2 ... (2n)^2} = \frac{(-1)^n}{2^{2n} n!^2} donc n=7n=7 convient.

Pour m0m \neq 0, on peut faire un raisonnement analogue (les calculs sont un peu plus compliqués).

On a ainsi trouvé une solution y 0y_0 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|x| \leq 1), on peut alors trouver toutes les solutions de l’équation différentielle (en posant y=y 0zy=y_0 z et en cherchant zz).

Exercice : faire de même pour les solutions de yxy=0y'{'}-xy=0 (fonctions de Airy).

4.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 +e ttdt f(x)=\int_x^{+\infty} \frac{e^{-t}}{t} \ dt On peut montrer que l’intégrale existe bien, car l’intégrand est positif et inférieur à e te^{-t} (qui admet e t-e^{-t} comme primitive, cette primitive ayant une limite en ++\infty). Pour trouver le développement asymptotique de ff en ++\infty, 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) = [e tt] x + x +e tt 2dt = e xx x +e tt 2dt = e xx([e tt 2] x + x +2e tt 3) = e xxe xx 2+ x +2e tt 3dt = ... = e x(1x1x 2+2x 3+...+(1) nn!x n+1) x +(1) n(n+1)!e tt n+2dt = S(x)+R(x) \begin{matrix} f(x)&=&[\frac{-e^{-t}}{t}]_x^{+\infty} - \int_x^{+\infty} \frac{-e^{-t}}{-t^2} \ dt \\ &=& \frac{e^{-x}}{x} - \int_x^{+\infty} \frac{e^{-t}}{t^2} \ dt \\ &=& \frac{e^{-x}}{x} - ([\frac{-e^{-t}}{t^2}]_x^{+\infty} - \int_x^{+\infty} \frac{-2e^{-t}}{-t^3}) \\ &=& \frac{e^{-x}}{x} - \frac{e^{-x}}{x^2} + \int_x^{+\infty} \frac{2e^{-t}}{t^3} \ dt \\ &=& ... \\ &=& e^{-x}\left(\frac{1}{x} - \frac{1}{x^2} + \frac{2}{x^3} + ... + \frac{(-1)^n n!}{x^{n+1}}\right) - \int_x^{+\infty} \frac{(-1)^n (n+1)!e^{-t}}{t^{n+2}} \ dt \\ &=& S(x) + R(x) \end{matrix} S(x)=e x(1x1x 2+2x 3+...+(1) nn!x n+1),R(x)= x +(1) n(n+1)!e tt n+2dt(4) S(x)=e^{-x} \left(\frac{1}{x} - \frac{1}{x^2} + \frac{2}{x^3} + ... + \frac{(-1)^n n!}{x^{n+1}}\right), \quad R(x)=- \int_x^{+\infty} \frac{(-1)^n (n+1)!e^{-t}}{t^{n+2}} \ dt \qquad (4) Le développement en séries est divergent puisque pour x>0x&gt;0 fixé et nn tendant vers l’infini lim n+n!x n+1=+ \lim_{n\rightarrow +\infty} \frac{n!}{x^{n+1}} = +\infty mais si xx est grand, au début la série semble converger, de manière très rapide : 1x>>1x 2>>2x 3 \frac{1}{x} &gt;&gt; \frac{1}{x^2} &gt;&gt; \frac{2}{x^3} On peut utiliser S(x)S(x) comme valeur approchée de f(x)f(x) pour xx grand si on sait majorer R(x)R(x) par un nombre suffisamment petit. On a |R(x)| x +(n+1)!e tx n+2=(n+1)!e xx n+2 | R(x) | \leq \int_x^{+\infty} \frac{(n+1)!e^{-t}}{x^{n+2}} = \frac{(n+1)!e^{-x}}{x^{n+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 e x/xe^{-x}/x. Pour xx fixé assez grand, il faut donc trouver un rang nn, s’il en existe un, tel que (n+1)!/x n+1<ε(n+1)!/x^{n+1}&lt;\epsilonε\epsilon est la précision relative que l’on s’est fixée. Par exemple, si x100x\geq 100, n=11n=11 convient pour ε=12!/100 12=5e16\epsilon=12!/100^{12}=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 xx).

Ce type de développement asymptotique peut être effectué pour d’autres fonctions du même type, par exemple x +e t 2dt, x +sin(t)tdt,... \int_x^{+\infty} e^{-t^2} \ dt, \quad \int_x^{+\infty} \frac{\sin(t)}{t} \ dt, \quad ...

Digression: calcul approché de la constante d’Euler γ\gamma
On peut montrer que lim n+u n,u n= k=1 n1kln(n)(5) \lim_{n\rightarrow +\infty} u_n, \quad u_n=\sum_{k=1}^{n}\frac{1}{k} - \ln(n) \qquad (5) existe (par exemple en cherchant un équivalent de u n+1u nu_{n+1}-u_n qui vaut 12n 2\frac{-1}{2n^2}) et on définit γ\gamma comme sa limite. Malheureusement, la convergence est très lente et cette définition n’est pas applicable pour obtenir la valeur de γ\gamma avec une très grande précision. Il y a un lien entre γ\gamma et la fonction exponentielle intégrale, plus précisément lorsque x0x\rightarrow 0, f(x)f(x) admet une singularité en ln(x)-\ln(x), plus précisément f(x)+ln(x)f(x)+\ln(x) admet un développement en séries (de rayon de convergence ++\infty), car : f(x)+ln(x) = x 1e t1tdt+ 1 +e ttdt = 0 1e t1tdt+ 1 +e ttdt 0 xe t1tdt \begin{matrix} f(x)+\ln(x)&=&\int_x^{1}\frac{e^{-t}-1}{t} \ dt + \int_1^{+\infty} \frac{e^{-t}}{t} \ dt \\ &=& \int_0^{1}\frac{e^{-t}-1}{t} \ dt + \int_1^{+\infty} \frac{e^{-t}}{t} \ dt - \int_0^{x} \frac{e^{-t}-1}{t} \ dt \end{matrix} Que vaut la constante du membre de droite : C= 0 1(e t1)1tdt+ 1 +e t1tdt C=\int_0^{1}(e^{-t}-1)\frac{1}{t} \ dt + \int_1^{+\infty} e^{-t} \frac{1}{t} \ dt Il se trouve que C=γC=-\gamma (voir plus bas une démonstration condensée) et donc : γ= 0 x1e ttdtf(x)ln(x)(6) \gamma= \int_0^{x} \frac{1-e^{-t}}{t} \ dt -f(x)-\ln(x) \qquad (6) Pour obtenir une valeur approchée de γ\gamma, il suffit donc de prendre un xx assez grand pour pouvoir calculer f(x)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=0x=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 nn suffisamment grand). Par exemple, on pose x=13x=13, on calcule f(13)f(13) par (4) avec n=13n=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/xn/x) et une erreur absolue inférieure à e 1313!/13 14=4e12e^{-13} 13!/13^{14}=4e-12

f(13)f(13) \approx exp(-13)*sum((-1)^n*n!/13.^(n+1),n=0..13)

puis on remplace dans (6), avec 0 x1e ttdt= n=0 (1) nx n+1(n+1)(n+1)! \int_0^{x} \frac{1-e^{-t}}{t} \ dt = \sum_{n=0}^{\infty} (-1)^n \frac{x^{n+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 R 50R_{50} 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 γ\gamma

-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 (5) qui définit γ\gamma 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)f(x) pour xx plus grand (déterminé par la précision qu’on peut obtenir par le développement aymptotique de ff).

On peut calculer π\pi 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 ff) et l’égalité (dont un schéma de preuve est aussi donné plus bas) 0 +sin(t)tdt=π2(7) \int_0^{+\infty} \frac{\sin(t)}{t} \ dt = \frac{\pi}{2} \qquad (7)

Calcul de CC (et preuve de (7)):
Pour cela on effectue une intégration par parties, cette fois en intégrant 1/t1/t et en dérivant l’exponentielle (moins 1 dans la première intégrale). C = 0 1(e t1)1tdt+ 1 +e t1tdt = [(e t1)ln(t)] 0 1+ 0 1ln(t)e tdt+[e tln(t)] 1 ++ 1 +ln(t)e tdt = 0 +ln(t)e tdt \begin{matrix} C&=&\int_0^{1}(e^{-t}-1)\frac{1}{t} \ dt + \int_1^{+\infty} e^{-t} \frac{1}{t} \ dt\\ &=&[(e^{-t}-1)\ln(t)]_0^1 +\int_0^1 \ln(t) e^{-t} \ dt + [e^{-t} \ln(t)]_1^{+\infty} +\int_1^{+\infty} \ln(t) e^{-t} \ dt \\ &=& \int_0^{+\infty} \ln(t) e^{-t} \ dt \end{matrix} 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 +t ne tdt n!= \int_0^{+\infty}t^n e^{-t} \ 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 tt et on dérive l’autre terme, puis t 2/2t^2/2, etc. C = [te tln(t)] 0 + 0 +te t(1tln(t))dt = 0 0 +e tdt+ 0 +te tln(t)dt = 1+[t 22e tln(t)] 0 + 0 +t 22e t(1tln(t))dt = 1 0 +t2e t+ 0 +t 22e tln(t)dt = 112+ 0 +t 22e tln(t)dt = ... = 112...1n+ 0 +t nn!e tln(t)dt = 112...1n+ln(n)+I n \begin{matrix} C&=&[te^{-t} \ln(t)]_0^{+\infty} - \int_0^{+\infty} t e^{-t}(\frac{1}{t}-\ln(t)) \ dt \\ &=& 0 - \int_0^{+\infty} e^{-t} \ dt + \int_0^{+\infty} t e^{-t} \ln(t) \ dt \\ &=& -1 + [\frac{t^2}{2}e^{-t} \ln(t)]_0^{+\infty} - \int_0^{+\infty} \frac{t^2}{2} e^{-t}(\frac{1}{t}-\ln(t)) \ dt \\ &=& -1 - \int_0^{+\infty} \frac{t}{2} e^{-t} + \int_0^{+\infty} \frac{t^2}{2} e^{-t} \ln(t) \ dt \\ &=& -1 - \frac{1}{2} + \int_0^{+\infty} \frac{t^2}{2} e^{-t} \ln(t) \ dt \\ &=& ...\\ &=& -1 - \frac{1}{2} - ... - \frac{1}{n} + \int_0^{+\infty} \frac{t^n}{n!} e^{-t} \ln(t) \ dt \\ &=& -1 - \frac{1}{2} - ... - \frac{1}{n} + \ln(n) + I_n \end{matrix} I n= 0 +t nn!e t(ln(t)ln(n))dt I_n=\int_0^{+\infty} \frac{t^n}{n!} e^{-t} (\ln(t)-\ln(n)) \ dt Pour déterminer I nI_n on fait le changement de variables t=nut=nu I n = 0 +(nu) nn!e nuln(u)ndu = n n+1n! 0 +e n(ln(u)u)ln(u)du \begin{matrix} I_n&=&\int_0^{+\infty} \frac{(nu)^n}{n!} e^{-nu} \ln(u) n\ du \\ &=& \frac{n^{n+1}}{n!} \int_0^{+\infty} e^{n(ln(u)-u)} \ln(u) \ du \end{matrix} Or en faisant le même changement de variables t=nut=nu : n!= 0 +t ne tdt=n n+1 0 +e n(ln(u)u)du n!= \int_0^{+\infty}t^n e^{-t} \ dt = n^{n+1} \int_0^{+\infty} e^{n(ln(u)-u)} \ du Donc I n= 0 +e n(ln(u)u)ln(u)du 0 +e n(ln(u)u)du I_n= \frac{\int_0^{+\infty} e^{n(ln(u)-u)} \ln(u) \ du} {\int_0^{+\infty} e^{n(ln(u)-u)} \ du} Lorsque nn tend vers l’infini, on peut montrer que I n0I_n \rightarrow 0, en effet les intégrales sont équivalentes à leur valeur sur un petit intervalle autour de u=1u=1, point où l’argument de l’exponentielle est maximal, et comme l’intégrand du numérateur a une amplitude ln(u)\ln(u) qui s’annule en u=1u=1, il devient négligeable devant le dénominateur. Finalement on a bien C=γC=-\gamma.

On peut remarquer qu’en faisant le même calcul que CC mais en remplacant e te^{-t} par e αte^{-\alpha t} pour (α)>0\Re(\alpha)&gt;0, donne limI n=ln(α)\lim I_n=-\ln(\alpha) (car le point critique où la dérivée de la phase s’annule est alors 1/α1/\alpha). Ceci peut aussi se vérifier pour α\alpha réel en faisant le changement de variables αt=u\alpha t=u 0 1(e αt1)1tdt+ 1 +e αt1tdt=γln(α) \int_0^{1}(e^{-\alpha t}-1)\frac{1}{t} \ dt + \int_1^{+\infty} e^{-\alpha t} \frac{1}{t} \ dt = -\gamma -\ln(\alpha) En faisant tendre α\alpha vers i-i, ln(α)-\ln(\alpha) tend vers ln(i)=iπ2\ln(i)=i\frac{\pi}{2} et on obtient 0 1(e it1)1tdt+ 1 +e it1tdt=γ+iπ2 \int_0^{1}(e^{it}-1)\frac{1}{t} \ dt + \int_1^{+\infty} e^{i t} \frac{1}{t} \ dt = -\gamma + i \frac{\pi}{2} dont la partie imaginaire nous donne (7), et la partie réelle une autre identité sur γ\gamma faisant intervenir la fonction cosinus intégral.

5  Polynômes : arithmétique, factorisation, interpolation

5.1  Arithmétique des polynomes: Bézout et applications

On considère les polynômes à une variable à coefficients dans \mathbb{R} ou \mathbb{C} ou \mathbb{Q}. Les algorithmes de base déjà évoqués sont l’évaluation en un point (méthode de Horner), l’addition, la soustraction, la multiplication et la division euclidienne de AA par B0B \neq 0 : A=BQ+R,deg(R)<deg(B) A = B Q + R, \quad \mbox{deg}(R)&lt; \mbox{deg}(B)

A l’aide de la division euclidienne, on peut calculer le PGCD de deux polynômes par l’algorithme d’Euclide. Nous allons présenter l’algorithme d’Euclide étendu (ou de Bézout)

Théorème 1   Étant donnés 2 polynômes AA et BB, il existe deux polynômes UU, VV tels que AU+BV=pgcd(A,B),deg(U)<deg(B),deg(V)<deg(A) AU + B V = \mbox{pgcd} (A,B) , \quad \mbox{deg}(U) &lt; \mbox{deg}(B), \mbox{deg}(V) &lt; \mbox{deg}(A)

Algorithme :
On construit en fait 3 suites (U n)(U_n), (V n)(V_n) et (R n)(R_n) telles que : AU n+BV n=R n AU_n + B V_n = R_n

Exemple :
A=x 31,B=x 2+1A=x^3-1, B=x^2+1, les rangs 0 et 1 sont donnés ci-dessus. Au rang 2, Q 0Q_0 est le quotient euclidien de AA par BB (fonction quo) donc xx, d’où U 2=1,V 2=x,R 2=x1 U_2=1, V_2= -x, R_2=-x-1 Puis on divise x 2+1x^2+1 par x1-x-1, quotient x+1-x+1, donc U 3=x1,V 3=1+x(x+1)=1+xx 2,R 3=2 U_3=x-1, V_3=1+x(-x+1)=1+x-x^2, R_3=2

Preuve de l’algorithme :
On montre facilement par récurrence que la relation AU n+BV n=R nAU_n + B V_n = R_n est conservée. Comme R nR_n est la suite des restes, le dernier reste non nul est bien le pgcd de AA et BB. D’autre part, examinons les degrés des V kV_k. Supposons que deg(A)(A) \geq deg(B)(B) (sinon on échange AA et BB). Au rang n=0n=0, V 0=0V_0=0 donc V 2=Q 0V 1V_2=-Q_0 V_1, aux rangs suivants le degré de Q nQ_n est non nul (car le degré de R n+1R_{n+1} est strictement inférieur au degré de R nR_n) On montre donc par récurrence que la suite des degrés de V nV_n est croissante et que : deg(V n+2)=deg(Q n)+deg(V n+1) \mbox{deg} (V_{n+2}) = \mbox{deg}(Q_n)+\mbox{deg}(V_{n+1}) Comme deg(Q n)(Q_n)=deg(R n)(R_n)-deg(R n+1)(R_{n+1}), on en déduit que deg(V n+2)+deg(R n+1)=deg(V n+1)+deg(R n)=...=deg(V 1)+deg(R 0)=deg(A) \mbox{deg}(V_{n+2})+\mbox{deg}(R_{n+1}) = \mbox{deg}(V_{n+1})+\mbox{deg}(R_{n}) = ... = \mbox{deg}(V_1)+\mbox{deg}(R_0 ) = \mbox{deg}(A) Donc si n+2n+2 est le rang du dernier reste non nul, V n+2=VV_{n+2}=V et degVV=degAA-degR n+1R_{n+1} est donc strictement inférieur au degré de AA (car R n+1R_{n+1}, l’avant-dernier reste non nul, est de degré plus grand ou égal à 1). On en déduit enfin que le degré de UU est strictement inférieur au degré de BB, car AU=RBVAU=R-BV, le degré de BVBV est strictement inférieur à celui de BB plus celui de AA.

L’identité de Bézout permet de résoudre plus générallement une équation du type Au+Bv=C Au+Bv=C A,B,CA,B,C sont trois polynômes donnés, à condition que CC soit divisible par le pgcd de AA et BB. L’ensemble des solutions s’obtient à partir d’une solution particulière U,VU,V de Bézout, notons c=C/gcd(A,B)c=C/\mbox{gcd}(A,B), on a alors A(cU)+B(cV)=cgcd(A,B)=C A(cU) + B(cV)=c\ \mbox{gcd}(A,B) = C et l’ensemble des solutions est donné par u=cUPB,v=cV+PAu=cU-PB, v=cV+PAPP est un polynôme quelconque. Si le degré de CC est plus petit que le degré de AA plus le degré de BB, il existe une solution “priviligiée”, on prend pour uu le reste de la division euclidienne de cUcU par BB, vv est alors le reste de la division euclidienne de cVcV par AA pour des raisons de degré.

Exemple : si on veut résoudre (x 31)u+(x 2+1)v=2x 2 (x^3-1)u+(x^2+1)v=2x^2 on multiplie U=x1U=x-1 et V=1+xx 2V=1+x-x^2 par x 2x^2 ce qui donne une solution u=x 2(x1),v=x 2(1+xx 2) u=x^2(x-1), \quad v=x^2(1+x-x^2) l’ensemble des solutions est de la forme u+P(x 2+1),vP(x 31) u+P(x^2+1), \quad v-P(x^3-1) et la solution priviligiée (de degrés minimaux) est x+1=rem(x 2(x1),x 2+1),x 2x+1=rem(x 2(1+xx 2),x 31) -x+1=\mbox{rem}(x^2(x-1),x^2+1), \quad x^2-x+1=\mbox{rem}(x^2(1+x-x^2),x^3-1)


L’identité de Bézout intervient dans de nombreux problèmes en particulier la décomposition en éléments simples d’une fraction rationnelle. Si le dénominateur DD d’une fraction se factorise en produit de 2 facteurs D=ABD=AB premiers entre eux, alors il existe deux polynômes uu et vv tels que N=Au+BvN=Au+Bv, donc ND=Au+BvAB=uB+vA \frac{N}{D} = \frac{Au+Bv} {AB} = \frac{u}{B} + \frac{v}{A} Si de plus N/DN/D est une fraction propre (degré de NN plus petit que celui de DD), alors u/Bu/B et v/Av/A sont encore des fractions propres (en calculant le reste de la division euclidienne pour uu et vv comme expliqué ci-dessus).

Par exemple : 2x 2(x 31)(x 2+1)=(x+1)(x 31)+(x 2x+1)(x 2+1)(x 31)(x 2+1)=x+1x 2+1+x 2x+1x 31 \frac{2x^2}{(x^3-1)(x^2+1)} = \frac{(-x+1)(x^3-1)+(x^2-x+1)(x^2+1)}{(x^3-1)(x^2+1)} = \frac{-x+1}{x^2+1} + \frac{x^2-x+1}{x^3-1}

Les applications sont diverses, citons

Il faut néanmoins savoir factoriser un polynôme, ce dont nous parlerons dans la section suivante.

Exercice : Calculer l’intégrale 1(x1)(x 2+1) \int \frac{1}{(x-1)(x^2+1)} en utilisant l’identité de Bézout pour décomposer la fraction rationnelle. Trouver à l’aide de cette décomposition le terme d’ordre nn du développement de Taylor de la fraction à intégrer, vérifier avec un logiciel de calcul formel que les termes d’ordre 0 à 3 sont corrects.


Una autre application est l’élimination dans les systèmes polynomiaux, par exemple considérons le système de 2 équations à 2 inconnues (intersection d’une ellipse et d’un cercle) : x 2+y 29=0,x 2+2y 22xy7=0 x^2+y^2-9=0, x^2+2y^2-2xy-7=0 En calculant les coefficients de Bézout des 2 polynômes en xx x 2+y 29x^2+y^2-9 et x 2+2y 22xy7x^2+2y^2-2xy-7 et en multipliant au besoin par le PPCM (plus grand commun multiple) des dénominateurs, on obtient à droite de l’équation de Bézout un polynôme ne dépendant que de yy et qui s’annule aussi aux solutions du système, on peut alors résoudre en yy (en factorisant) puis en xx. Ici par exemple ce polynome est 5y 432y 2+45y^4-32y^2+4.



Cette méthode se systématise, le polynome obtenu par élimination d’une variable est appelé résultant.



5.2  Factorisation des polynômes

Soit PP un polynôme de degré non nul. Factoriser PP n’a pas une signification unique, tout dépend d’une part si on veut une factorisation exacte ou approchée, et d’autre part quels seront les types des coefficients de la factorisation (complexes, réels, entiers).

5.2.1  Multiplicité des racines.

On dit que rr est une racine de multiplicité kk de PP si P(x)=(xr) kQP(x)=(x-r)^k Q et Q(r)0Q(r)\neq 0.

En faisant le développement de Taylor de PP en rr à l’ordre degré de PP, on voit que cela équivaut à : P(r)=P(r)=...=P [k1](r)=0,P [k](r)0 P(r)=P'(r)=...=P^{[k-1]}(r)=0, \quad P^{[k]}(r) \neq 0 En particulier si P(r)=0P(r)=0, on peut factoriser PP par XrX-r.

On peut donc détecter les racines de multiplicité supérieure à 1 en cherchant un facteur commun à PP et PP', en effet xrx-r divisera PP et PP'.

Théorème 2   Si PP et PP' sont premiers entre eux (pgcd = 1), alors les racines de PP sont simples (de multiplicité 1).

P:=x^4-1; gcd(P,P’)

On a aussi le résultat suivant :

Proposition 3   Soit PP un polynôme à coefficients réels ou complexes, et soit Q=P/Q=P/pgcd(P,P)(P,P'). Alors PP et QQ ont les mêmes racines et les racines de QQ sont simples.

Preuve : si P(z)=C i(zz i) m iP(z)=C\prod_i (z-z_i)^{m_i} avec z iz_i les racines distinctes de PP, alors pgcd(P,P)= i(zz i) m i1(P,P')=\prod_i (z-z_i)^{m_i-1} et Q= i(zz i)Q=\prod_i (z-z_i).


Ce résultat est très utile si le polynôme donné est á coefficients exacts, car les méthodes numériques d’approximation de racines ne fonctionnent bien que pour des polynômes à racines simples.

Il existe un algorithme (dû à Yun) qui permet d’écrire un polynome quelconque comme produit de polynômes dont les racines sont simples en effectuant uniquement des calculs de PGCD de polynomes.

yun(P):= { 
  local W,Y,G,res; 
  W:=P;  
  Y:=diff(W,x);  
  res:=NULL;  
  while(true){ 
    if (Y==0) { 
      return res[1..size(res)-1],W;  
    };  
  G:=gcd(Y,W);  
  res:=res,G;  
  W:=normal(W/G);  
  Y:=normal(Y/G-diff(W,x));  
  };  
}:;

onload


L’instruction sqrfree ou équivalente de votre logiciel de calcul formel effectue cette décomposition.

5.2.2  Factorisation dans \mathbb{C}.

Reste maintenant à trouver des racines! On a le :

Théorème 4   (d’Alembert)
Soit
PP un polynome de degré non nul, alors PP admet au moins une racine complexe.

On peut alors factoriser PP par XrX-r si rr est la racine, et recommencer avec le quotient, d’où le corollaire.

Théorème 5   Soit PP un polynome de degré nn non nul, alors PP admet nn racines complexes (comptées avec multiplicité) x 1,...,x nx_1,...,x_n, on a donc : P(X)=a n j=1 n(Xx j)(8) P(X)=a_n \prod_{j=1}^n (X-x_j) \qquad (8) a na_n est le coefficient dominant de PP.

Démonstration du théorème de d’Alembert :
On va montrer que le minimum de la valeur absolue de PP est atteint en un nombre complexe puis que ce minimum est forcément nul. Soit P(x)=a nx n+...+a 0,a n0 P(x)=a_n x^n + ... + a_0, \quad a_n \neq 0 Lorsque |x||x| tend vers l’infini, |P(x)||P(x)| tend vers l’infini, en effet P(x)=a nx n(1+a n1a n1x+...+a 0a n1x n) |x|a nx n P(x)= a_n x^n ( 1 + \frac{a_{n-1}}{a_n} \frac{1}{x} + ... + \frac{a_0}{a_n} \frac{1}{x^n}) \approx_{|x|\rightarrow \infty} a_n x^n plus précisément il existe R>0R&gt;0 tel que si |x|>R|x|&gt;R alors |P(x)|>|a n||x| n/2|P(x)|&gt;|a_n| |x|^n/2. Quitte à augmenter RR on peut donc supposer que |P(x)|>|P(0)||P(x)|&gt;|P(0)| si |x|>R|x|&gt;R, donc il existe un complexe x 0x_0 qui réalise le minimum de |P||P| sur \mathbb{C} (ce minimum est en fait le minimum pour |x|R|x|\leq R). On va montrer par l’absurde que ce minimum est nul (donc que x 0x_0 est la racine cherchée). Supposons donc que P(x 0)0P(x_0) \neq 0. On fait le développement de Taylor de PP en x 0x_0 à l’ordre nn=degré de PP, donc le développement n’a pas de reste : P(x)P(x 0)=(xx 0)P(x 0)+..+(xx 0) nP [n](x 0)n! P(x) - P(x_0) = (x-x_0) P'(x_0) + .. + (x-x_0)^n \frac{P^{[n]}(x_0)}{n!} Comme PP n’est pas constant, l’un des termes du membre de droite est non nul, soit kk l’indice du premier terme non nul, on a alors : P(x)=P(x 0)+(xx 0) kP [k](x 0)k!+o((xx 0) k) P(x) = P(x_0) + (x-x_0)^k \frac{P^{[k]}(x_0)}{k!} + o((x-x_0)^k) Comme P(x 0)0P(x_0) \neq 0, on peut le factoriser en : P(x)=P(x 0)(1+(xx 0) kP [k](x 0)P(x 0)k!+o((xx 0) k) P(x) = P(x_0)( 1 + (x-x_0)^k \frac{P^{[k]}(x_0)}{P(x_0) k!} + o((x-x_0)^k) on pose alors x=x 0+twx=x_0+t www est une racine kk-ième (cela existe dans \mathbb{C}) de (P [k](x 0)P(x 0)k!) 1 \left(\frac{-P^{[k]}(x_0)}{P(x_0) k!} \right)^{-1} on a alors : P(x)=P(x 0)(1t k+o(t k+1)) P(x)=P(x_0)(1-t^k+o(t^{k+1})) lorsque tt est positif, suffisamment petit, on a 0<1t k+o(t k+1<10 &lt; 1-t^k+o(t^{k+1} &lt; 1, donc |P(x)|<|P(x 0)||P(x)|&lt;|P(x_0)|, ce qui est absurde (x 0x_0 réalisant le minimum de PP sur \mathbb{C}).

Remarque :
Si on développe la relation (8), on obtient des relations entre les coefficients du polynome et les racines, par exemple : a n1=a nj=1 n(x j),...,a 0=a n j=1 n(x j), a_{n-1}=a_n \sum{j=1}^n (-x_j), ...,\quad a_0 = a_n \prod_{j=1}^n (-x_j),

5.2.3  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 dans la section suivante qui permet de connaitre le signe de PP'{'} 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 QQ de PP par XrX-r (par l’algorithme de Horner), puis on calcule les racines du quotient QQ (qui sont des racines de PP).

Un problème pratique apparait alors, c’est que rr n’est pas exact donc le quotient QQ non plus, au fur et à mesure du calcul des racines de PP, on perd de plus en plus de précision. Il existe une amélioration simple, si rr' est une racine approchée de QQ, alors elle est racine approchée de PP et on a toutes les chances qu’elle soit suffisamment proche d’une racine de PP pour que le théorème s’applique, on effectue alors 1 ou 2 itérations de Newton avec rr' mais pour PP (et non QQ) afin d’améliorer sa précision comme racine de PP. Si PP est à coefficients exacts, pour éviter les erreurs de calculs en approché, on remplace rr par un rationnel proche et on effectue ces itérations de Newton en calcul exact (entre deux itérations, pour des raisons d’efficacité, on peut remplacer rr par un rationnel proche à la précision attendue à cette étape mais dont le numérateur et dénominateur ont une écriture moins longue).

5.2.4  Localisation d’une racine complexe près d’une racine approchée

On peut montrer à postériori des estimations sur la distance entre une racine approchée et la racine la plus proche d’un polynôme, plus précisément cette distance est inférieure ou égale au degré du polynôme multiplié par le module de P/PP/P' en la racine approchée.

Proposition 6   Le disque de centre zz et de rayon degre(P)|P(z)/P(z)|(P) |P(z)/P'(z)| contient au moins une racine de PP.

En effet PP(z)= k1zz k \frac{P'}{P}(z)= \sum_{k} \frac{1}{z-z_k} donc si |zz k|>|z-z_k| &gt; degre(P)|P(z)/P(z)|(P) |P(z)/P'(z)| pour toutes les racines z kz_k alors |PP(z)|>|P(z)/P(z)||\frac{P'}{P}(z)| &gt; |P'(z)/P(z)| contradiction.

En combinant des méthodes d’approximation numérique des racines d’un polynôme et ce résultat on peut ainsi localiser très précisément des racines complexes de polynôme. C’est un cas où une méthode hybride exacte-approchée donne des résultats excellents!

Remarque : pour des polynômes de degré nn dépassant la dizaine, on utilise des méthodes permettant de trouver toutes les racines à la fois (diagonalisation numérique de la matrice companion) plutôt que la méthode de Newton qui les détermine une par une. Lorsque le degré dépasse quelques centaines, il devient nécessaire d’utiliser des approximations multi-précision ou des résutats de fonction à variable complexes.

5.2.5  Localisation des racines réelles : Sturm

Pour factoriser un polynôme à coefficients réels, on commence par le factoriser dans \mathbb{C}. On observe ensuite que si rr est une racine complexe non réelle de PP, alors son conjugué l’est aussi (il suffit de prendre le conjugue de la relation P(r)=0P(r)=0) et avec la même multiplicité (les dérivées successives de PP étant aussi à coefficients réels). On regroupe alors les facteurs correspondant à des racines complexes conjuguées : (Xr)(Xr¯)=X 2(r+r¯)X+rr¯=X 22(r)X+|r| 2 (X-r)(X-\overline{r}) = X^2 - (r+\overline{r})X+r\overline{r} = X^2 - 2 \Re(r) X + |r|^2 Finalement, on a le :

Théorème 7   La factorisation d’un polynôme à coefficients réels sur \mathbb{R} donne un produit de facteurs de degré 1 (correspondant à des racines réelles) et de degré 2 (correspondant à des paires de racines complexes conjuguées)

Il existe un algorithme utilisant l’algorithme de calcul du PGCD de PP et PP' qui permet de déterminer le nombre de racines réelles d’un polynôme PP sans racine multiple sur \mathbb{R} ou dans un intervalle de RR.

Théorème 8   On définit la suite de polynômes A 0=P,A 1=P,...,A k,0A_0=P, A_1=P', ..., A_k,0 en prenant l’opposé du reste de la division euclidienne des deux précdents : A i=A i+1Q i+2A i+2(9) A_{i} = A_{i+1} Q_{i+2} - A_{i+2} \qquad (9) Soit A kA_k, le dernier reste non nul, c’est un polynôme constant puisque PP n’a pas de racine multiple. On définit s(a)s(a) comme étant le nombre de changements de signes de la suite A i(a)A_i(a) en ignorant les 0. Alors le nombre de racines réelles de A 0=PA_0=P sur l’intervalle ]a,b]]a,b] est égal à s(a)s(b)s(a)-s(b).

Exemple :
Quel est le nombre de racines réelles de P=x 3+x+1P=x^3+x+1 sur [2,2][-2,2]? sur [0,2][0,2]?
On a donc A 0=x 3+x+1,A 1=P=3x 2+1,A 2=rem(A 0,A 1,x)=23x1,A 3=314 A_0=x^3+x+1, \quad A_1=P'=3x^2+1, \quad A_2=-\mbox{rem}(A_0,A_1,x)=-\frac{2}{3} x-1, \quad A_3=-\frac{31}{4} En x=2x=-2 on obtient la suite 9,13,1/3,31/4-9,13,1/3,-31/4 (2 changements de signe), en x=2x=2 on obtient la suite 11,13,7/3,31/411,13,-7/3,-31/4 (1 changement de signe), il y a donc 1 racine réelle entre -2 et 2. En x=0x=0 on obtient la suite 1,1,1,31/41,1,-1,-31/4 (1 changement de signe) donc la racine réelle est entre -2 et 0.


Preuve
On considére la suite des signes en un point : elle ne peut contenir deux 0 successifs (sinon toute la suite vaudrait 0 en ce point en appliquant (9), or A kA_k est constant non nul). Elle ne peut pas non plus contenir ...,+,0,+,... ni ...,-,0,-,... à cause de la convention de signe sur les restes de (9). Donc si bb est une racine de A iA_i pour 0<i<k0&lt;i&lt;k, alors en bb on a soit ...,-,0,+,... soit ...,+,0,-,... . Regardons le premier cas (le deuxième cas se traite de manière analogue), pour xx proche de bb, on va avoir ...,-,-,+,... ou ...,-,+,+,... dans les 2 cas la contribution au nombre de changements de signe est constant (égal à 1).

Comme A kA_k est constant, seules les racines de A 0=PA_0=P sont susceptibles de faire varier ss. Comme A 1=PA_1=P', le sens de variations de A 0A_0 au voisinage d’une racine de A 0A_0 est déterminé par le signe de A 1A_1, donc lorsque xx augmente en traversant une racine rr de PP, il y a deux possibilités soit PP est croissant et on passe de -,+,... à +,+,..., soit PP est décroissant et on passe +,-,... à -,-,.... Dans les deux cas, on diminue ss d’une unité.

Application :
Si il n’existe pas de racines réelles dans un intervalle donné, alors le polynôme garde un signe constant sur cet intervalle, que l’on peut déterminer en calculant la valeur du polynôme en un point de cet intervalle. On peut ainsi établir dans certains cas que la méthode de Newton pour trouver une racine d’un polynôme convergera.

Par exemple pour le polynôme P=3x 510x 3+30x 2x45P=3x^5-10x^3+30x^2-x-45, on a P=60(x 3x+1)P'{'}=60(x^3-x+1), est positif sur +\mathbb{R}^+ (exercice : calculer la suite de Sturm correspondante pour le vérifier). On vérifie que P(1)<0P(1)&lt;0 et P(1)>0P'(1)&gt;0 donc il existe une racine r>1r&gt;1 telle que P(r)>0P'(r)&gt;0, toute valeur de départ de Newton supérieure à rr assure la convergence.

Remarque :
On peut aussi déterminer les racines réelles d’un polynôme à coefficients rationnels en faisant uniquement des calculs exacts par dichotomie. Cette méthode de localisation des racines réelles se généralise d’ailleurs au cas complexe. On peut ainsi déterminer les racines complexes d’un polynôme à coefficients complexes rationnels de manière déterministe à la précision voulue (cf. Eisermann).

5.2.6  Localisation des racines réelles : règle des signes de Descartes

Proposition 9   Soit P=a nx n+...+a 1x+a 0[x]P=a_nx^n+...+a_1x+a_0 \in \mathbb{R}[x], r(P)r(P) le nombre de racines de PP sur *+\mathbb{R}^{*+}, et σ(P)\sigma(P) le nombre de changements de signe de la suite des coefficients a n,...,a 0a_n,...,a_0 de PP. Alors r(P)σ(P)r(P) \leq \sigma(P).

La preuve se fait par récurrence. Pour n=0n=0 on a bien 000 \leq 0. Pour n>0n&gt;0 on applique l’hypothèse de récurrence à PP' et on utilise le fait que 1+le nombre de racines de PP' sur un intervalle ]a,b[]a,b[ est supérieur ou égal au nombre de racines de PP sur [a,b][a,b]. Il y a deux cas possibles σ(P)=σ(P)\sigma(P')=\sigma(P) et σ(P)=σ(P)1\sigma(P')=\sigma(P)-1. Dans le deuxième cas, comme r(P)+1r(P)r(P')+1\geq r(P), alors r(P)r(P)+1σ(P)+1=σ(P)r(P) \leq r(P')+1 \leq \sigma(P')+1=\sigma(P). Dans le premier cas, il faut trouver une racine supplémentaire pour PP'. Pour cela, on regarde ce qui se passe en 0 +0^+. Supposons que a 0>0a_0&gt;0 pour fixer les idées, comme a 0a_0 est du même signe que a 1a_1 (ou que a 2a_2 si a 1=0a_1=0 etc.), PP' est positif en 0 +0^+ donc PP croit en 0 +0^+, donc doit atteindre un maximum local avant la première racine de PP, ce maximum local est une racine de PP'.

On peut même montrer que r(P)r(P) et σ(P)\sigma(P) sont de même parité. En particulier si σ(P)=1\sigma(P)=1, on a r(P)=1r(P)=1 car le coefficient dominant et de plus bas degré non nul de PP sont de signes contraires, donc la valeur en 0 et la limite en ++\infty aussi.

On déduit de la règle des signes de Descartes un critère indiquant si un polynôme possède 0, 1 ou plus de racines sur ]0,1[]0,1[ en se ramenant par changement de variables x1/xx \leftarrow 1/x à ]1,+[]1,+\infty[ puis par xx+1x \leftarrow x+1 à ]0,+[]0,+\infty[.

On prend alors le polynôme de départ et une majoration sur les racines, par exemple C=|P| /|lcoeff(P)|C=|P|_{\infty}/|\mbox{lcoeff}(P)|. Si zz est une racine positive de PP, alors z/Cz/C est une racine positive de Q(x):=P(Cx)Q(x):=P(Cx), on est donc ramené a chercher des intervalles d’isolation de racines de QQ dans ]0,1[]0,1[. S’il y a 0 ou 1 changement de signe, on conclut. Sinon, on teste si 1/21/2 est racine puis on cherche dans ]0,1/2[]0,1/2[ et dans ]1/2,1[]1/2,1[, ce qui revient au problème précédent sur ]0,1[]0,1[ avec les polynômes R(x)=Q(x/2)R(x)=Q(x/2) et S(x)=Q((x+1)/2)S(x)=Q((x+1)/2).

On fait de même sur P(x=x)P(x=-x) pour les racines réelles négatives.

fonction descartes(P)
  local l,res,j;
  l:=coeffs(P);
  l:=remove(0,l);
  res:=0;
  pour j de 0 jusque dim(l)-2 faire
    si l[j]*l[j+1]<0 alors res++; fsi;
  fpour;
  return res;
ffonction:;
  
fonction isole01(P,a,b)
  // renvoie une liste d'intervalles d'isolation
  local n,m,Q,R,res;
  global x;
  Q:=x^degree(P)*P(x=1/x);
  R:=Q(x=x+1);
  n:=descartes(R);
  si n=0 alors return NULL; fsi;
  si n=1 alors return [a,b]; fsi;
  Q:=numer(P(x=x/2));
  R:=numer(P(x=(x+1)/2));
  m:=(a+b)/2;
  si P(x=1/2)=0 alors res:=(a+b)/2; sinon res:=NULL; fsi;
  res:=res,isole01(Q,a,m),isole01(R,m,b)
  retourne res;
ffonction:;

fonction isole(P) // racines positives de P
  local l,M;
  l:=coeffs(P);
  M:=ceil(maxnorm(l)/abs(l[0]))+1;
  P:=P(x=M*x);
  return isole01(P,0,M);
ffonction:;

onload
Racines de PP dans *+\mathbb{R}^{*+}

il y a un nombre pair de racines positives. Racines de PP dans ]0,1[]0,1[

Donc une racine dans ]0,1[]0,1[.

Isolation des racines positives de PP: on commence par se ramener à chercher les racines de QQ sur ]0,1[]0,1[ (il faudra les multiplier par MM à la fin)



Il faut découper en deux :

donc il y a une racine de QQ dans ]0,1/2[]0,1/2[, donc de PP dans ]0,5/2[]0,5/2[

donc il y a une racine de QQ dans ]1/2,1[]1/2,1[, donc de PP dans ]5/2,5[]5/2,5[.

QQ n’est pas nul en 1/21/2, on conclut donc qu’il existe 2 racines positives pour PP. Vérification avec isole :

5.2.7  Factorisation exacte

Soit PP un polynôme à coefficients entiers. Lorsqu’on demande à un logiciel de calcul formel de factoriser PP, par défaut il ne calcule pas les racines complexes approchées, mais renvoie une factorisation exacte, sous forme de produit de facteurs à coefficients entiers. Les degrés des facteurs peuvent être plus grand que 2. Par exemple x 4+x+1x^4+x+1 ne peut pas être factorisé en produit de polynômes à coefficients entiers (bien qu’il ait 2 facteurs de degré 2 dans \mathbb{R} et 4 de degré 1 dans \mathbb{C}).

Commencons par une méthode simple de calcul des racines rationnelles de PP (les racines rationnelles correspondent à des facteurs entiers de degré 1 de la forme qXpqX-p de PP). Soit x=p/qx=p/q une racine rationnelle écrite sous forme de fraction irréductible de P=a nX n+...+a 0P=a_n X^n+...+a_0, on a alors 0=P(pq)=a np nq n+a n1p n1q n1+...+a 0=a np n+a n1p n1q+...+a 1pq n1+a 0q nq n 0 = P(\frac{p}{q}) = a_n \frac{p^n}{q^n} + a_{n-1} \frac{p^{n-1}}{q^{n-1}} + ... + a_0 = \frac{a_n p^n + a_{n-1} p^{n-1}q+...+a_1 p q^{n-1}+ a_0 q^n}{q^n} Donc : p(a np n1+a n1p n1q+...+a 1q n1)=a 0q n p(a_n p^{n-1} + a_{n-1} p^{n-1}q+...+a_1 q^{n-1}) = - a_0 q^n et pp divise donc a 0q na_0 q^n. Comme p/qp/q est irréductible, cela entraine que pp divise a 0a_0. De même qq divise a na_n. Il suffit donc de tester quelles sont les racines de PP parmi toutes les fractions irréductibles de la forme un diviseur de a 0a_0 sur un diviseur de a na_n (attention à ne pas oublier les diviseurs négatifs!).

Exemple: racines rationnelles de 2x 2+3x+1=02x^2+3x+1=0. On a pp divise 11 donc vaut 1 ou -1, qq divise 22 donc vaut 1 ou 2. On teste donc 1, -1, 1/2, -1/2. On obtient ici la factorisation complète du polynome (les racines sont -1 et -1/2) 2x 2+3x+1=2(x+1)(x+1/2) 2x^2+3x+1=2(x+1)(x+1/2)

Remarques :

Pour déterminer les facteurs à coefficients entiers de plus grand degré, il n’existe pas de méthode aussi simple. On peut calculer des valeurs approchées des racines complexes et essayer de créer des paquets de racines complexes, puis tester si a n rpaquet(Xr)a_n\prod_{r \in \mbox{paquet}} (X-r) est à coefficients entier (aux erreurs d’arrondi près). Par exemple si on calcule les racines complexes approchées de x 6+2x 3x 2+1x^6+2x^3-x^2+1, on pourra composer un facteur de degré 3 à coefficients entiers en rassemblant les racines de x 3+x+1x^3+x+1. Les logiciels de calcul formel utilisent des algorithmes modulaires et pp-adiques (consistant à factoriser le polynome modulo pp).

5.3  Approximation polynomiale

É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 x 0,...,x nx_0,...,x_n et à trouver une majoration de la différence entre la fonction et le polynôme. Le polynome interpolateur de Lagrange répond à cette question.

Soit donc x 0,...,x nx_0,...,x_n des réels distincts et y 0,...,y ny_0,...,y_n les valeurs de la fonction à approcher en ces points (on posera y j=f(x j)y_j=f(x_j) pour approcher la fonction ff). On cherche donc PP tel que P(x j)=y iP(x_j)=y_i pour j[0,n]j \in [0,n].

Commencons par voir s’il y a beaucoup de solutions. Soit PP et QQ deux solutions distinctes du problème, alors PQP-Q est non nul et va s’annuler en x 0,...,x nx_0, ...,x_n donc possède n+1n+1 racines donc est de degré n+1n+1 au moins. Réciproquement, si on ajoute à PP un multiple du polynome A= j=0 n(Xx j)A=\prod_{j=0}^n (X-x_j), 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+1n+1 multiple de AA.

Nous allons maintenant construire une solution particulière de degré au plus nn. Si n=0n=0, on prend P=x 0P=x_0 constant. On procède ensuite par récurrence. Pour construire le polynôme correspondant à x 0,...,x n+1x_0,...,x_{n+1} on part du polynoôme P nP_n correspondant à x 0,...,x nx_0,...,x_{n} et on lui ajoute un multiple réel de AA P n+1=P n+α n+1 j=0 n(Xx j) P_{n+1}=P_n+ \alpha_{n+1} \prod_{j=0}^n (X-x_j) Ainsi on a toujours P n+1(x j)=y jP_{n+1}(x_j)=y_j pour j=0,..nj=0,..n, on calcule maintenant α n+1 \alpha_{n+1} pour que P n+1(x n+1)=y n+1P_{n+1}(x_{n+1})=y_{n+1}. En remplacant avec l’expression de P n+1P_{n+1} ci-dessus, on obtient P n(x n+1)+α n+1 j=0 n(x n+1x j)=y n+1 P_n(x_{n+1})+ \alpha_{n+1} \prod_{j=0}^n (x_{n+1}-x_j) = y_{n+1} Comme tous les x jx_j sont distincts, il existe une solution unique : α n+1=y n+1P n(x n+1) j=0 n(x n+1x j) \alpha_{n+1}=\frac{y_{n+1}-P_n(x_{n+1})}{\prod_{j=0}^n (x_{n+1}-x_j)}

On a donc prouvé le :

Théorème 10   Soit n+1n+1 réels distincts x 0,...,x nx_0,...,x_n et n+1n+1 réels quelconques y 0,...,y ny_0,...,y_n. Il existe un unique polynôme PP de degré inférieur ou égal à nn, appelé polynome de Lagrange, tel que : P(x i)=y i P(x_i)=y_i

Exemple : déterminons le polynome de degré inférieur ou égal à 2 tel que P(0)=1,P(1)=2,P(2)=1P(0)=1, P(1)=2, P(2)=1. On commence par P 0=1P_0=1. Puis on pose P 1=P 0+α 1X=1+α 1XP_1=P_0+ \alpha_{1}X=1+ \alpha_{1}X. Comme P(1)=2=1+α 1P(1)=2=1+ \alpha_{1} on en tire α 1=1 \alpha_{1}=1 donc P 1=1+XP_1=1+X. Puis on pose P 2=P 1+α 2X(X1)P_2=P_1+ \alpha_{2}X(X-1), on a P 2(2)=3+2α 2=1P_2(2)=3+2 \alpha_{2}=1 donc α 2=1 \alpha_{2}=-1, finalement P 2=1+XX(X1)P_2=1+X-X(X-1).

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

Théorème 11   Soit ff une fonction n+1n+1 fois dérivable sur un intervalle I=[a,b]I=[a,b] de \mathbb{R}, x 0,...,x nx_0,...,x_n des réels distincts de II. Soit PP le polynome de Lagrange donné par les x jx_j et y j=f(x j)y_j=f(x_j). Pour tout réel xIx \in I, il existe un réel ξ x[a,b]\xi_x \in [a,b] (qui dépend de xx) tel que : f(x)P(x)=f [n+1](ξ x)(n+1)! j=0 n(xx j)(10) f(x)-P(x) = \frac{f^{[n+1]}(\xi_x)}{(n+1)!} \prod_{j=0}^n(x-x_j) \qquad (10)

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

Preuve du théorème : Si xx est l’un des x jx_j l’égalité est vraie. Soit C=(f(x)P(x))/ j=0 n(xx j) C=(f(x)-P(x))/\prod_{j=0}^n(x-x_j) on considère maintenant la fonction : g(t)=f(t)P(t)C j=0 n(tx j) g(t)=f(t)-P(t) - C \prod_{j=0}^n(t-x_j) elle s’annule en x jx_j pour jj variant de 0 à nn ainsi qu’en xx suite au choix de la constante CC, donc gg s’annule au moins n+2n+2 fois sur l’intervalle contenant les x jx_j et xx, donc gg' s’annule au moins n+1n+1 fois sur ce même intervalle, donc gg'{'} s’annule au moins nn fois, etc. et finalement g [n+1]g^{[n+1]} s’annule une fois au moins sur cet intervalle. Or g [n+1]=f [n+1]C(n+1)! g^{[n+1]} = f^{[n+1]} - C (n+1)! car PP est de degré inférieur ou égal à nn et j=0 n(xx j)x n+1 \prod_{j=0}^n(x-x_j) - x^{n+1} est de degré inférieur ou égal à nn. Donc il existe bien un réel ξ x\xi_x dans l’intervalle contenant les x jx_j et xx tel que C=f [n+1](ξ x)(n+1)! C=\frac{f^{[n+1]}(\xi_x)}{(n+1)!}

Illustration : approcher la fonction f(x)=x*sin(x)cos(x)f(x)=x*\sin(x)-\cos(x) sur [1,2][-1,2] avec un pas de 0.5.





Graphiquement, on ne distingue pas f(x)f(x) et PP sur l’intervalle d’interpolation. Tracons la différence :

Dès qu’on sort de l’intervalle, l’erreur explose. Comparons avec la majoration d’erreur sur l’intervalle



La dérivée 7-ième est donc majorée par 10 sur [1,2][-1,2]. On trace la majoration de l’erreur avec l’erreur.

On observe qu’avec des points d’interpolation équirépartis, l’erreur est plus grande aux bords de l’intervalle qu’au milieu. Pour rendre la majoration de l’erreur uniforme, il faudrait utiliser plus de points vers les bords de l’intervalle, c’est ce qu’on fait sur [1,1][-1,1] en utilisant les racines du polynôme de Tchebycheff de 1ère espèce T n+1T_{n+1} défini par T k(cos(x))=cos(kx)T_k(\cos(x))=\cos(kx) Ce polynôme de degré kk est uniformément bornée par 1 sur [1,1][-1,1]. Et si un polynôme UU a le même degré et coefficient dominant que T kT_k alors sa norme est au moins 1, car si elle était strictement plus petite que 1, la différence T kUT_k-U serait du signe de T kT_k lorsque T k=±1T_k=\pm 1 et changerait donc de signe kk fois sur l’intervalle [1,1][-1,1] ce qui est impossible pour un polynôme de degré au plus k1k-1.

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(xx 0)+...+α n(xx 0)...(xx n1) = α 0+(xx 0)(α 1+(xx 1)(α 2+...+(xx n2)(α n1+(xx n1)α n)...)) \begin{matrix} P(x) &=& \alpha_0 + \alpha_1 (x-x_0) + ... + \alpha_n (x-x_0)...(x-x_{n-1}) \\ &=& \alpha_0 + (x-x_0)( \alpha_1 + (x-x_1)(\alpha_2 + ... + (x-x_{n-2})(\alpha_{n-1}+(x-x_{n-1}) \alpha_n)...)) \end{matrix} ce qui permet de le calculer rapidement une fois les α i\alpha_i connus. On observe que α 0=f(x 0),α 1=f(x 1)f(x 0)x 1x 0 \alpha_0=f(x_0), \quad \alpha_1=\frac{f(x_1)-f(x_0)}{x_1-x_0} On va voir que les α k\alpha_k peuvent aussi se mettre sous forme d’une différence. On définit les différences divisées d’ordre nn par récurrence f[x i]=f(x i),f[x i,...,x k+i+1]=f[x i+1,...,x k+i+1]f[x i,...,x k+i]x k+i+1x i f[x_i]=f(x_i), \quad f[x_i,...,x_{k+i+1}]= \frac{f[x_{i+1},...,x_{k+i+1}]-f[x_i,...,x_{k+i}]}{x_{k+i+1}-x_i} On va montrer que α k=f[x 0,...,x k]\alpha_k=f[x_0,...,x_k]. C’est vrai au rang 0, il suffit donc de le montrer au rang k+1k+1 en l’admettant au rang kk. Pour cela on observe qu’on peut construire le polynôme d’interpolation en x 0,...,x k+1x_0,...,x_{k+1} à partir des polynômes d’interpolation P kP_k en x 0,...,x kx_0,...,x_k et Q kQ_k en x 1,...,x k+1x_1,...,x_{k+1} par la formule : P k+1(x)=(x k+1x)P k+(xx 0)Q kx k+1x 0 P_{k+1}(x)= \frac{(x_{k+1}-x)P_k + (x-x_0)Q_k}{x_{k+1}-x_0} en effet on vérifie que P k+1(x i)=f(x i)P_{k+1}(x_i)=f(x_i) pour i[1,k]i\in [1,k] car P k(x i)=f(x i)=Q k(x i)P_k(x_i)=f(x_i)=Q_k(x_i), et pour i=0i=0 et i=k+1i=k+1, on a aussi P k+1(x 0)=f(x 0)P_{k+1}(x_0)=f(x_0) et P k+1(x k+1)=f(x k+1)P_{k+1}(x_{k+1})=f(x_{k+1}). Or α k+1\alpha_{k+1} est le coefficient dominant de P k+1P_{k+1} donc c’est la différence du coefficient dominant de Q kQ_k et de P kP_k divisée par x k+1x 0x_{k+1}-x_0, c’est-à-dire la définition de f[x 0,...,x k+1]f[x_0,...,x_{k+1}] en fonction de f[x 1,...,x k+1]f[x_1,...,x_{k+1}] et f[x 0,...,x k]f[x_0,...,x_{k}].

Exemple : on reprend P(0)=1,P(1)=2,P(2)=1P(0)=1, P(1)=2, P(2)=1. On a x i f[x i] f[x i,x i+1] f[x 0,x 1,x 2] 0 framebox1 (21)/(10)=framebox1 1 2 (11)/(20)=framebox1 (12)/(21)=1 2 1 \begin{array}{cccc} x_i & f[x_i] & f[x_i,x_{i+1}] & f[x_0,x_1,x_2] \\ 0 & \framebox{1} & & \\ & & (2-1)/(1-0)=\framebox{1} & \\ 1 & 2 & & (-1-1)/(2-0)=\framebox{-1} \\ & & (1-2)/(2-1)=-1 & \\ 2 & 1 & & \\ \end{array} donc P(x)=framebox1+(x0)(framebox1+(x1)(framebox1))=1+x(2x)P(x)=\framebox{1}+(x-0)(\framebox{1}+(x-1)(\framebox{-1}))=1+x(2-x).

On peut naturellement utiliser l’ordre que l’on souhaite pour les x ix_i, en observant que le coefficient dominant de PP ne dépend pas de cet ordre, on en déduit que f[x 0,...,x k]f[x_0,...,x_k] est indépendant de l’ordre des x ix_i, on peut donc à partir du tableau ci-dessus écrire PP par exemple avec l’ordre 2,1,0, sous la forme P(x)=1+(x2)(1+(x1)(1))=1+(x2)(x) P(x)=1+(x-2)(-1+(x-1)(-1))=1+(x-2)(-x)

6  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+(n1)h,a+nh=b[a,b]=[a,a+h] + [a+h,a+2h]+...[a+(n-1)h,a+nh=b, où h=(ba)/nh=(b-a)/n est le pas de la subdivision, et sur chaque subdivision, on approche l’aire sous la courbe.

6.1  Les rectangles et les trapèzes

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

Exemple : calcul de la valeur approchée de 0 1t 3dt\int_0^1 t^3 dt (on en connait la valeur exacte 1/4=0.251/4=0.25) par ces méthodes en subdivisant [0,1][0,1] en 10 subdivisions (pas h=1/10h=1/10), donc α=j/10\alpha=j/10 et β=(j+1)/10\beta=(j+1)/10 pour jj variant de 0 à 9. Pour les rectangles à gauche, on obtient sur une subdivision f(α)=(j/10) 3f(\alpha)=(j/10)^3 que l’on multiplie par la longueur de la subdivision soit h=1/10h=1/10 : 110 j=0 9(j10) 3=81400=0.2025 \frac{1}{10} \sum_{j=0}^9 (\frac{j}{10})^3 = \frac{81}{400} = 0.2025 Pour les rectangles à droite, on obtient 110 j=1 10(j10) 3=121400=0.3025 \frac{1}{10} \sum_{j=1}^{10} (\frac{j}{10})^3 = \frac{121}{400} = 0.3025 Pour le point milieu f((α+β)/2)=f((j/10+(j+1)/10)/2)=f(j/10+1/20)f((\alpha+\beta)/2)=f((j/10+(j+1)/10)/2)=f(j/10+1/20) 110 j=0 9(j10+120) 3=199/800=0.24875 \frac{1}{10} \sum_{j=0}^9 (\frac{j}{10}+\frac{1}{20})^3 = 199/800 = 0.24875 Enfin pour les trapèzes, l’aire du trapèze délimité par l’axe des xx, les verticales y=αy=\alpha, y=βy=\beta et les points sur ces verticales d’ordonnées respectives f(α)f(\alpha) et f(β)f(\beta) vaut hf(α)+f(β)2 h \frac{f(\alpha)+f(\beta)}{2} donc 110 j=0 9((j10) 3+(j+110) 3)=101400=0.2525 \frac{1}{10} \sum_{j=0}^9 \left( (\frac{j}{10})^3 +(\frac{j+1}{10})^3 \right) = \frac{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érallement, les formules sont donc les suivantes : rectangle gauche = h j=0 n1f(a+jh) rectangle droit = h j=1 nf(a+jh) point milieu = h j=0 n1f(a+jh+h2) trapezes = h(f(a)+f(b)2+ j=1 n1f(a+jh))(11) \begin{matrix} \mbox{rectangle gauche} & = & h \sum_{j=0}^{n-1} f(a+jh) \\ \mbox{rectangle droit} & = & h \sum_{j=1}^{n} f(a+jh) \\ \mbox{point milieu} & = & h \sum_{j=0}^{n-1} f(a+jh+\frac{h}{2}) \\ \mbox{ trapezes } & = & h \left(\frac{f(a)+f(b)}{2}+\sum_{j=1}^{n-1} f(a+jh) \right) \end{matrix} \qquad (11) h=(ba)/nh=(b-a)/n est le pas de la subdivision, nn 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 ff est continument dérivable, de dérivée majorée par une constante M 1M_1 sur [a,b][a,b], en faisant un développement de Taylor de ff en α\alpha, on obtient | α βf(t)dt α βf(α)dt|=| α βf(θ t)(tα)dt|M 1 α β(tα)dt=M 1(βα) 22 |\int_{\alpha}^{\beta} f(t) dt - \int_{\alpha}^{\beta} f(\alpha) dt | = | \int_{\alpha}^{\beta} f'(\theta_t)(t-\alpha) dt | \leq M_1 \int_{\alpha}^{\beta} (t-\alpha) dt = M_1\frac{(\beta-\alpha)^2}{2} Ainsi dans l’exemple, on a M 1=3M_1=3, l’erreur est donc majorée par 0.0150.015 sur une subdivision, donc par 0.150.15 sur les 10 subdivisions.

Pour le point milieu, on fait le développement en (α+β)/2(\alpha+\beta)/2 à l’ordre 2, en supposant que ff est deux fois continument dérivable : | α βf(t) α βf(α+β2)| = | α βf(α+β2)(tα+β2)dt + α βf(θ t)2(tα+β2) 2| M 222 α+β2 β(tα+β2) 2dt M 2(βα) 324 \begin{matrix} |\int_{\alpha}^{\beta} f(t) - \int_{\alpha}^{\beta} f(\frac{\alpha+\beta}{2}) | &= &| \int_{\alpha}^{\beta} f'(\frac{\alpha+\beta}{2})(t-\frac{\alpha+\beta}{2}) dt \\ & & + \int_{\alpha}^{\beta} \frac{f'{'}(\theta_t)}{2}(t-\frac{\alpha+\beta}{2})^2 | \\ &\leq & \frac{M_2}{2} 2 \int_{\frac{\alpha+\beta}{2}}^{\beta} (t-\frac{\alpha+\beta}{2})^2 dt \\ & \leq & M_2\frac{(\beta-\alpha)^3}{24} \end{matrix} Dans l’exemple, on a M 2=6M_2=6, donc l’erreur sur une subdivision est majorée par 0.25e30.25e-3, donc sur 10 subdivisions par 0.25e2=0.00250.25e-2=0.0025.

Pour les trapèzes, la fonction gg dont le graphe est le segment reliant [α,f(α)][\alpha,f(\alpha)] à [β,f(β)][\beta,f(\beta)] est f(α)+(tα)/(βα)f(β)f(\alpha)+(t-\alpha)/(\beta-\alpha)f(\beta), c’est en fait un polynome de Lagrange, si ff est deux fois continument dérivable, on peut donc majorer la différence entre ff et gg en utilisant (10), on intègre la valeur absolue ce qui donne | α βf(t)dt α βg(t)dt| α β|f(ξ x)2(xα)(xβ)|M 2(βα) 312 |\int_{\alpha}^{\beta} f(t) dt - \int_{\alpha}^{\beta} g(t) dt | \leq \int_{\alpha}^{\beta} |\frac{f'{'}(\xi_x)}{2} (x-\alpha)(x-\beta)| \leq M_2 \frac{(\beta-\alpha)^3}{12} M 2M_2 est un majorant de |f||f'{'}| sur [a,b][a,b].

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

Lorsque hh 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 hh va être grande, plus la convergence sera rapide lorsque hh sera petit, avec toutefois une contrainte fixée par la valeur de M kM_k, borne sur la dérivée kk-ième de ff (plus kk est grand, plus M kM_k est grand en général). Nous allons voir dans la suite comment se comporte cette puissance de hh en fonction de la facon dont on approche ff.

6.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(t)dtI(f)= j=1 kw jf(y j) \int_{\alpha}^{\beta} f(t) dt \approx I(f)=\sum_{j=1}^k w_j f(y_j) où les y jy_j sont dans l’intervalle [α,β][\alpha,\beta], par exemple équirépartis sur [α,β][\alpha,\beta]. On utilise aussi la définition : α βf(t)dtI(f)=(βα) j=1 kw˜ jf(y j) \int_{\alpha}^{\beta} f(t) dt \approx I(f)= (\beta-\alpha)\sum_{j=1}^k \tilde{w}_j f(y_j) On prend toujours jw j=βα\sum_j w_j=\beta-\alpha (ou jw˜ j=1\sum_j \tilde{w}_j=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 nn si il y a égalité ci-dessus pour tous les polynômes de degré inférieur ou égal à nn et non égalité pour un polynôme de degré n+1n+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érallement, si on approche ff par son polynôme d’interpolation de Lagrange en n+1n+1 points (donc par un polynôme de degré inférieur ou égal à nn), on obtient une méthode d’intégration d’ordre au moins nn.

Si une méthode est d’ordre nn avec des w j0w_j\geq 0 et si ff est n+1n+1 fois continument dérivable, alors sur une subdivision, on a : | α βfI(f)|M n+1(βα) n+2(n+1)!(1n+2+1)(12) |\int_{\alpha}^{\beta} f-I(f)| \leq M_{n+1} \frac{(\beta-\alpha)^{n+2}}{(n+1)!} (\frac{1}{n+2}+1) \qquad (12)

En effet, on fait le développement de Taylor de ff par exemple en α\alpha à l’ordre nn f(t) = T n(f)+(tα) n+1(n+1)!f [n+1](θ t), T n(f) = f(α)+(tα)f(α)+...+(tα) nn!f [n](α) \begin{matrix} f(t)&=&T_{n}(f)+\frac{(t-\alpha)^{n+1}}{(n+1)!} f^{[n+1]}(\theta_t),\\ T_{n}(f)&=&f(\alpha)+(t-\alpha)f'(\alpha)+...+ \frac{(t-\alpha)^{n}}{n!} f^{[n]}(\alpha) \end{matrix} Donc | α βf α βT n(f)| α β(tα) n+1(n+1)!|f [n+1](θ t)|[M n+1(tα) n+2(n+2)!] α β |\int_{\alpha}^{\beta} f- \int_{\alpha}^{\beta} T_{n}(f)| \leq \int_{\alpha}^{\beta} \frac{(t-\alpha)^{n+1}}{(n+1)!} |f^{[n+1]}(\theta_t)| \leq \left[ M_{n+1} \frac{(t-\alpha)^{n+2}}{(n+2)!} \right]_\alpha^\beta De plus, |I(f)I(T n(f))|=|I(f [n+1](θ t)(tα) n+1(n+1)!)| j=1 k|w j|M n+1(y jα) n+1(n+1)! j=1 k|w j|M n+1(βα) n+1(n+1)! \begin{matrix} |I(f) -I(T_n(f))| =|I\left( f^{[n+1]}(\theta_t) \frac{(t-\alpha)^{n+1}}{(n+1)!} \right)| & \leq & \sum_{j=1}^k |w_j| M_{n+1} \frac{(y_j-\alpha)^{n+1}}{(n+1)!} \\ & \leq & \sum_{j=1}^k |w_j| M_{n+1} \frac{(\beta-\alpha)^{n+1}}{(n+1)!} \end{matrix} Donc comme la méthode est exacte pour T n(f)T_n(f), on en déduit que | α βfI(f)| = | α βf α βT n(f)+I(T n(f))I(f)| | α βf α βT n(f)|+|I(T n(f))I(f)| M n+1(βα) n+2(n+2)!+ j=1 k|w j|M n+1(βα) n+1(n+1)! \begin{matrix} |\int_{\alpha}^{\beta} f-I(f)| &= &|\int_{\alpha}^{\beta} f-\int_{\alpha}^{\beta} T_n(f)+I(T_n(f))- I(f)| \\ &\leq& |\int_{\alpha}^{\beta} f-\int_{\alpha}^{\beta} T_n(f)|+|I(T_n(f))- I(f)|\\ &\leq & M_{n+1} \frac{(\beta-\alpha)^{n+2}}{(n+2)!} + \sum_{j=1}^k |w_j| M_{n+1} \frac{(\beta-\alpha)^{n+1}}{(n+1)!} \end{matrix} Si les w j0w_j\geq 0, alors j=1 k|w j|= j=1 kw j=βα\sum_{j=1}^k |w_j|=\sum_{j=1}^k w_j=\beta-\alpha et on obtient finalement (12)

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

Théorème 1   Pour une méthode d’ordre nn à coefficients positifs et une fonction ff n+1n+1 fois continument dérivable  | a bfI(f)|M n+1h n+12 n+1(n+1)!(ba)(1(n+2)+1)|\int_{a}^{b} f-I(f)| \leq M_{n+1} \frac{h^{n+1}}{2^{n+1}(n+1)!} (b-a) (\frac{1}{(n+2)}+1)

On observe que cette majoration a la bonne puissance de hh 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 nn.

6.3  Simpson

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

factor(int(lagrange([a,(a+b)/2,b],[fa,fm,fb]),x=a..b))

qui donne la formule sur une subdivision I(f)=h6(f(α)+4f(α+β2)+f(β)) I(f) = \frac{h}{6} (f(\alpha)+4f(\frac{\alpha+\beta}{2}) + f(\beta)) et sur [a,b][a,b] : I(f)=h6(f(a)+f(b)+4 j=0 n1f(a+jh+h2)+2 j=1 n1f(a+jh))(13) I(f) = \frac{h}{6} \left( f(a)+f(b)+ 4 \sum_{j=0}^{n-1} f(a+jh+\frac{h}{2}) + 2 \sum_{j=1}^{n-1} f(a+jh) \right) \qquad (13) Si on intègre t 3t^3 sur [0,1][0,1] en 1 subdivision par cette méthode, on obtient 16(0+412 3+1)=14 \frac{1}{6} (0+ 4 \frac{1}{2^3} + 1)=\frac{1}{4} c’est-à-dire le résultat exact, ceci est aussi vérifié pour ff polynome de degré inférieur ou égal à 2 puisque l’approximation de Lagrange de ff est alors égale à ff. 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 t 4t^4 sur [0,1][0,1] n’est pas exacte). On peut même améliorer (cf. par exemple Demailly) la constante générale de la section précédente pour la majoration de l’erreur en : | a bfI(f)|h 42880(ba)M 4 |\int_a^b f - I(f)| \leq \frac{h^4}{2880} (b-a) M_4 Cette méthode nécessite 2n+12n+1 évaluations de ff (le calcul de ff est un point étant presque toujours l’opération la plus couteuse en temps d’une méthode de quadrature), au lieu de nn pour les rectangles et le point milieu et n+1n+1 pour les trapèzes. Mais on a une majoration en h 4h^4 au lieu de h 2h^2 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 ff n’a pas la régularité suffisante (ou si M 4M_4 est trop grand).

6.4  Newton-Cotes

On peut généraliser l’idée précédente, découper la subdivision [α,β][\alpha,\beta] en nn parts égales et utiliser le polynôme d’interpolation en ces n+1n+1 points x 0=α,x 1,...,x n=βx_0=\alpha, x_1, ..., x_n=\beta. Ce sont les méthodes de Newton-Cotes, qui sont d’ordre nn 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)=(βα) j=0 nw˜ jf(x j) I(f)=(\beta-\alpha)\sum_{j=0}^n \tilde{w}_j f(x_j) De plus les w˜ j\tilde{w}_j sont universels (ils ne dépendent pas de la subdivision), parce qu’on peut faire le changement de variables x=α+t(βα)x=\alpha+t(\beta-\alpha) dans l’intégrale et le polynôme d’interpolation et donc se ramener à [0,1][0,1].

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

int(lagrange(seq(j/4,j,0,4),[f0,f1,f2,f3,f4]),x=0..1)

et on lit les coefficients de f0 à f4 qui sont les w˜ 0\tilde{w}_0 à w˜ 4\tilde{w}_4: 7/90, 32/90, 12/90, 32/90, 7/90. 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 | a bfI(f)|M 62 66!(1+17)h 6(ba) |\int_a^b f -I(f)| \leq \frac{M_6}{2^6 6!}(1+\frac{1}{7}) h^6 (b-a) elle peut être améliorée pour cette méthode précise en | a bfI(f)|M 61935360h 6(ba) |\int_a^b f -I(f)| \leq \frac{M_6}{1935360} h^6 (b-a)

En pratique, on ne les utilise pas très souvent, car d’une part pour n8n\geq 8, les w jw_j ne sont pas tous positifs, et d’autre part, parce que la constante M nM_n 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) 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][a,b], puis 2, puis 2 22^2, ..., et on élimine les puissances de hh du reste fI(f)\int f-I(f) en utilisant un théorème d’Euler-Mac Laurin qui montre que le développement asymptotique de l’erreur en fonction de hh ne contient que des puissances paires de hh). De plus, on peut être amené à faire varier le pas hh en fonction de la plus ou moins grande régularité de la fonction.

6.5  En résumé

Intégration sur [a,b][a,b], hh pas d’une subdivision, M kM_k majorant de la dérivée kk-ième de la fonction sur [a,b][a,b]

 formuleLagrange degréordreerreur
rectangles(11)00M 1h(ba)/2M_1 h (b-a)/2
point milieu(11)01M 2h 2(ba)/24M_2 h^2 (b-a)/24
trapèzes(11)11M 2h 2(ba)/12M_2 h^2 (b-a)/12
Simpson(13)23M 4h 4(ba)/2880M_4 h^4 (b-a)/2880

7  Algèbre linéaire

7.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 AA et BB reliés par une relation Ax=BAx=B, cette relation reste alors vraie au cours et donc après la réduction.

7.1.1  L’algorithme

L’algorithme est le suivant:

  1. on initialise c=1c=1 et l=1l=1, cc désigne le numéro de colonne cc à réduire, et ll le numéro de ligne à partir duquel on cherche un “pivot” (au début ll et cc valent donc 1, en général les 2 augmentent de 1 à chaque itération)
  2. Si cc ou ll est plus grand que le nombre de colonnes ou de lignes on arrête.
  3. Si la colonne cc n’a que des coefficients nuls à partir de la ligne ll, on incrémente le numéro de colonne cc 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 ll. Puis on effectue pour toutes les lignes sauf ll ou pour toutes les lignes à partir de l+1l+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 L jL ja jca lcL l L_j \leftarrow L_j -\frac{a_{jc}}{a_{lc}} L_l On incrémente cc et ll de 1 et on passe à l’étape 2.

7.1.2  Efficacité de l’algorithme

Si la matrice possède LL lignes et CC colonnes, le nombre maximal d’opérations pour réduire une ligne est CC divisions, CC multiplications, CC soustractions, donc 3C3C opérations arithmétiques de base. Il y a L1L-1 lignes à réduire à chaque étape et min(L,C)(L,C) étapes à effectuer, on en déduit que le nombre maximal d’opérations pour réduire une matrice est 3LC3LCmin(L,C)(L,C). Pour une matrice carrée de taille nn, cela fait 3n 33n^3 opérations.

7.1.3  Erreurs d’arrondis du pivot de Gauss

Comme |a jc||a lc||a_{jc}| \leq |a_{lc}|, 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 2 n2^n (où nn 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=52n=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 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, cf. l’exemple suivant.

Exemple
Soit à résoudre le système linéaire εx+1.0y=1.0,x+2.0y=3.0 \epsilon x + 1.0 y = 1.0 , \quad x + 2.0 y = 3.0 avec ε=2 54\epsilon =2^{-54} (pour une machine utilisant des doubles pour les calculs en flottant, plus générallement on choisira ε\epsilon tel que (1.0+3ε)1.0(1.0+3\epsilon)-1.0 soit indistinguable de 0.0).
Si on résoud le système exactement, on obtient x=1/(12ε)x=1/(1-2\epsilon) (environ 1) et y=(13ε)/(12ε)y=(1-3\epsilon)/(1-2\epsilon) (environ 1). Supposons que l’on n’utilise pas la stratégie du pivot partiel, on prend alors comme pivot ε\epsilon, donc on effectue la manipulation de ligne L 2L 21/εL 1L_2 \leftarrow L_2 - 1/\epsilon L_1 ce qui donne comme 2ème équation (2.01.0/ε)y=3.01.0/ε(2.0-1.0/\epsilon)y=3.0-1.0/\epsilon. 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/ε(-1.0/\epsilon)y=-1.0/\epsilon d’où y=1.0y=1.0, qui sera remplacé dans la 1ère équation, donnant εx=1.01.0y=0.0\epsilon x = 1.0-1.0y=0.0 donc x=0.0x=0.0.
Inversement, si on utilise la stratégie du pivot partiel, alors on doit échanger les 2 équations L 2=L 1L_2'=L_1 et L 1=L 2L_1'=L_2 puis on effectue L 2L 2εL 1L_2 \leftarrow L_2' - \epsilon L_1', ce qui donne (1.02.0ε)y=1.03.0ε(1.0-2.0\epsilon) y = 1.0 - 3.0 \epsilon , remplacée en raison des erreurs d’arrondi par 1.0*y=1.01.0*y=1.0 donc y=1.0y=1.0, puis on remplace yy dans L 1L_1' ce qui donne x=3.02.0y=1.0x=3.0-2.0y=1.0.
On observe dans les deux cas que la valeur de yy est proche de la valeur exacte, mais la valeur de xx 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. Le traitement du sujet à ce niveau est un peu difficile, cela fait intervenir le nombre de conditionnement de la matrice AA 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).

7.2  Applications de Gauss

7.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)(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)(1,0,-1), (0,1,2) (remarque: une réduction sous-diagonale aurait suffi).

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

7.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=bAx=bxx est la solution de l’équation, et à la fin de la réduction A=IA=I.

Par exemple pour résoudre le système {cccx+2y+3z = 5 4x+5y+6z = 0 7x+8y = 1 \left\{ \begin{matrix}{ccc} x + 2y + 3z &=& 5 \\ 4x + 5y + 6z &=& 0 \\ 7x + 8y &=& 1 \end{matrix} \right. on réduit sous forme échelonnée la matrice [[1,2,3,5],[4,5,6,0],[7,8,0,1]], 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/3x=-9, y=8, z=-2/3.

7.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 Ax=BA x=Bxx est l’inverse de la matrice de départ, et en fin de réduction A=IA=I.

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

7.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)(2,-1,0,0), la 4ème colonne donne le deuxième vecteur (2,0,2,1)(-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).

7.2.6  La méthode de factorisation LULU

Nous ne la développons pas à ce niveau, elle permet d’écrire une matrice AA comme produit de deux matrices triangulaire inférieures et supérieures, ce qui permet de ramener la résolution de système à la résolution de deux systèmes triangulaires.

7.3  Réduction exacte des endomorphismes

On calcule le polynome caractéristique ou le polynome minimal, on le factorise, et on calcule ensuite le noyau de AλIA-\lambda I pour les λ\lambda racines. Il existe des méthodes évitant le calcul de noyau, méthode de Fadeev-Laguerre-Souriau que nous ne présentons pas ici.

7.3.1  Polynome caractéristique

On peut le calculer en développant le déterminant detAλIA-\lambda I, mais il est plus efficace de le calculer par interpolation. Soit AA une matrice carrée de taille nn, on sait que son polynome caractéristique est un polynome de degré nn, il suffit de connaitre sa valeur en n+1n+1 points distincts, on calcule donc n+1n+1 déterminants detAλIA-\lambda I en remplaçant λ\lambda par sa valeur (il y a plus de déterminants à calculer mais ce sont des déterminants sans paramètre λ\lambda donc beaucoup plus simple à calculer), ce qui permet de reconstruire le polynome caractéristique par interpolation de Lagrange.

Exercice : pour [[1,-1],[2,4]], calculer det(AλI(A-\lambda I) en λ=0,1,2\lambda=0,1,2 puis le polynome d’interpolation, vérifier que c’est bien le polynome caractéristique.

Il faut effectuer n+1n+1 calculs de déterminants, ce qui nécessite O(n 4)O(n^4) opérations. Il existe des méthodes plus efficaces, par exemple le calcul du polynome minimal probabiliste présenté plus bas (O(n 3)O(n^3) opérations).

7.3.2  Polynome minimal

Définition 1   Le polynome minimal d’une matrice AA est un polynôme MM de degré minimal tel que M(A)=0M(A)=0 et de coefficient dominant égal à 1. Un tel polynome divise tous les polynomes tels que P(A)=0P(A)=0, il divise le polynome caractéristique de AA et il a les mêmes racines que le polynome caractéristique.

Preuve:
D’abord MM divise tous les polynomes tels que P(A)=0P(A)=0, car si RR désigne le reste de la division de PP par MM alors R(A)=(PQM)(A)=P(A)Q(A)M(A)=0R(A)=(P-QM)(A)=P(A)-Q(A)M(A)=0, donc RR est nul car son degré est plus petit que celui de MM.

En particulier le polynome minimal divise le polynome caractéristique CC, car C(A)=0C(A)=0 (on peut montrer que C(A)=0C(A)=0 en faisant le produit de la matrice AλIA-\lambda I par sa comatrice, on obtient le déterminant fois l’identité, soit C(λ)IC(\lambda)I. Comme C(λ)IC(A)C(\lambda)I-C(A) peut se factoriser par λIA\lambda I -A en appliquant (14) à chaque monome de CC, on en déduit que C(A)C(A) se factorise par λIA\lambda I -A, donc C(A)=0C(A)=0 en regardant les termes de plus haut degré de ces polynomes en λ\lambda à coefficients matriciels).

Montrons enfin que les racines du polynome caractéristique sont racines du polynome minimal. En effet soit λ\lambda une racine du polynome caractéristique alors AλIA-\lambda I n’est pas inversible. Or M(A)M(λ)IM(A)-M(\lambda)I se factorise par AλIA-\lambda I car A kλ kI=(AλI) j=0 k1λ k1jA j(14) A^k-\lambda^k I=(A-\lambda I) \sum_{j=0}^{k-1} \lambda^{k-1-j} A^j \qquad (14) donc M(A)M(λ)IM(A)-M(\lambda)I ne peut pas être inversible. Comme M(A)=0M(A)=0 on en déduit que M(λ)IM(\lambda)I n’est pas inversible donc M(λ)=0M(\lambda)=0, λ\lambda est une racine de MM. Donc si le polynome caractéristique n’a pas de racines multiples, il est égal au polynome minimal.

Pour calculer MM, on peut chercher une relation de degré minimal entre les puissances de AA, en les voyant comme des vecteurs à n 2n^2 composantes (ce qui revient à aplatir en un long vecteur tous les coefficients de la matrice). Cela revient à calculer le noyau de l’application linéaire dont les colonnes sont les coefficients des puissances de AA (de 0 à nn), en gardant le premier vecteur obtenu par l’algorithme calcul du noyau ci-dessus.

Cette méthode est toutefois couteuse, car il faut réduire une matrice ayant n 2n^2 lignes et n+1n+1 colonnes. Il existe une autre méthode moins couteuse et qui marche presque toujours. Elle consiste à calculer le polynome minimal de AA par rapport à un vecteur vv c’est-à-dire le polynome de degré minimal (et coefficient dominant 1) tel que M v(A)v=0M_v(A)v=0. Comme M(A)=0M(A)=0, on a M(A)v=0M(A)v=0, donc M vM_v divise MM qui divise le polynome caractéristique. Si par chance, on trouve que M vM_v est de degré nn, alors M vM_v sera égal à MM et au polynome caractéristique. On fait donc le calcul du noyau de l’application linéaire dont les colonnes sont les A jvA^jv pour jj variant de 0 à nn. Si l’on trouve un espace de dimension 1, alors M vM_v est de degré nn et on a simultanément le polynome minimal et caractéristique avec le polynome correspondant à ce vecteur du noyau. Si le degré n’est pas nn, on peut essayer un ou quelques autres vecteurs, et faire le PPCM des polynomes minimaux obtenus. Si on obtient un polynome de degré nn on conclut, sinon on peut tester si ce polynome évalué en AA est nul, ce sera alors le polynome minimal.

Exemple, on reprend la matrice [[1,-1],[2,4]], et comme vecteur aléatoire v=(1,0)v=(1,0), on a Av=(1,1)Av=(1,-1) et A(Av)=(1,5)A(Av)=(-1,-5). On calcule donc le noyau de la matrice [[1,1,-1],[0,-1,-5]] (on écrit en colonnes v,Av,A 2vv, Av, A^2v), on trouve que (6,5,1)(-6,5,-1) engendre le noyau, donc le polynome minimal relatif au vecteur vv est (au signe près) 6+5xx 2-6+5x-x^2. Comme il est de degré maximal 2, c’est le polynome minimal et caractéristique.

7.4  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λIA-\lambda I), lorsqu’on veut calculer quelques valeurs propres on préfère utiliser des méthodes itératives directement sur AA ce qui évite la propagation des erreurs d’arrondi.

7.4.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 AA soient x 1,...,x nx_1,...,x_n avec |x 1||x 2|...|x n1|<|x n||x_1| \leq |x_2| \leq ... \leq |x_{n-1}| &lt; |x_n| et soient e 1,...,e ne_1,...,e_n une base de vecteurs propres correspondants. On choisit un vecteur aléatoire vv et on calcule la suite v n=Av n1=A nvv_n=Av_{n-1}=A^n v . Si vv a pour coordonnées V 1,...,V n)V_1,...,V_n) dans la base propre, alors v n= j=1 nV jx j ne j=x n nw n,w n=V j(x jx n) ne j v_n = \sum_{j=1}^n V_j x_j^n e_j = x_n^n w_n, \quad w_n=\sum V_j \left(\frac{x_j}{x_n}\right)^n e_j L’hypothèse que x nx_n est l’unique valeur propre de module maximal entraine alors que lim n+w n=V ne n\lim_{n \rightarrow +\infty} w_n = V_n e_n puisque la suite géométrique de raison x j/x nx_j/x_n converge vers 0. Autrement dit, si V n0V_n\neq 0 (ce qui a une probabilité 1 d’être vrai pour un vecteur aléatoire), v nv_n est équivalent à V nx n ne nV_n x_n^n e_n. Lorsque nn est grand, v nv_n est presque colinéaire au vecteur propre e ne_n (que l’on peut prendre comme v nv_n divisé par sa norme), ce que l’on détecte en testant si v n+1v_{n+1} et v nv_n sont presques colinéaires, et de plus le facteur de colinéarité entre v n+1v_{n+1} et v nv_n est presque x nx_n, la valeur propre de module maximal.

Exercice : tester la convergence de A nvA^nv vers l’espace propre associé à λ=3\lambda=3 pour la matrice [[1,-1],[2,4]] et le vecteur (1,0)(1,0).

Lorsqu’on applique cette méthode a une matrice réelle, il peut arriver qu’il y ait deux valeurs propres conjuguées de module maximal. Le même type de raisonnement montre que pour nn grand, v n+2v_{n+2} est presque colinéaire à l’espace engendré par v nv_n et v n+1v_{n+1}, la relation v n+2+xv n+1+x 2v n=0 v_{n+2}+ x v_{n+1} + x^2 v_n=0 permet de calculer les valeurs propres.

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

7.4.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 xx, on peut appliquer la méthode précédente à la matrice (AxI) 1(A-xI)^{-1}. En effet, les valeurs propres de cette matrice sont les (x ix) 1(x_i-x)^{-1} dont la norme est maximale lorsqu’on se rapproche de x ix_i.

7.4.3  Elimination des valeurs propres trouvées

Si la matrice AA est symétrique, et si e ne_n est un vecteur propre normé écrit en colonne, on peut considérer la matrice B=Ax ne ne n tB=A-x_n e_n e_n^t qui possède les mêmes valeurs propres et mêmes vecteurs propres que AA avec même multiplicité, sauf x nx_n qui est remplacé par 0. En effet les espaces propres de AA sont orthogonaux entre eux, donc Be n=x ne nx ne ne n te n=0,Be k=x ke kx ne ne n te k=x ke k Be_n=x_ne_n -x_n e_n e_n^t e_n = 0, Be_k = x_k e_k - x_n e_n e_n^t e_k = x_k e_k On peut donc calculer la 2ème valeur propre (en valeur absolue), l’éliminer et ainsi de suite.

Si la matrice AA n’est pas symétrique, on peut utiliser une technique analogue lorsque 0 n’est pas valeur propre de AA (on peut s’y ramener en ajoutant à AA un multiple de l’identité). En effet on peut construire un vecteur propre de BB pour une valeur propre x k0x_k \neq 0 à partir d’un vecteur propre de BB, en cherchant yy tel que tel que B(e kye n)=x k(e kye n) B(e_k-ye_n)=x_k(e_k-ye_n) On obtient pour le membre de gauche : Be kyBe n=Be k=(Ax ne ne n t)e k=x ke kx ne n.e ke n Be_k - yB e_n =Be_k= (A-x_n e_n e_n^t)e_k = x_k e_k - x_n e_n.e_k e_n et pour le membre de droite x ke kyx ke n x_k e_k - y x_k e_n d’où l’équation yx k=x ne n.e k y x_k = x_n e_n.e_k Néanmoins cette méthode n’est pas stable, en particulier si la valeur propre e ke_k est proche de 0, car les vecteurs propres se rapprochent alors tous de e ne_n.

8  Guide rapide KhiCAS sur calculatrices

8.1  Installation

8.1.1  Installation Casio Graph 90

(Ceci s’applique dans une certaine mesure aux Casio Graph 35eii). Rendez-vous ici pour installer. Documentation Casio

8.1.2  Installation TI Nspire CX/CX II

Vous pouvez installer KhiCAS sur TI Nspire CX/CX II à condition de ne pas avoir installé la toute dernière mise à jour de l’OS de TI. Rendez-vous ici pour installer. Documentation Nspire.

8.1.3  Installation Numworks N0110 non verrouillée.

Pour les Numworks, il faut un modèle N0110 (juillet 2019 ou plus tard) et ne surtout pas avoir fait de mise à jour vers Epsilon 16 ou supérieur. Rendez-vous sur cette page avec un navigateur compatible web-USB. Documentation Numworks.

8.2  Le shell KhiCAS

Si KhiCAS n’est pas lancé :

Pour éteindre la calculatrice, il faut quitter KhiCAS sur les TI Nspire (taper doc deux fois avant ctrl on) et les Casio (taper MENU avant shift AC/ON).

Attention, sur les TI Nspire et Numworks, le shell de KhiCAS peut être configuré pour utiliser un interpréteur MicroPython, le bandeau du bas apparait alors en jaune, et sur les TI Nspire un interpréteur Javascript, le bandeau du bas apparait alors en orange. Pour changer d’interpréteur, sur les Nspire utiliser la touche calculatrice (entre esc et tab) puis 2. Sur les Numworks shift ) 8.

Vous êtes maintenant dans le shell de KhiCAS, vous pouvez y taper des commandes Xcas en vous aidant des menus qui apparaissent dans le bandeau et de raccourcis claviers (F1 à F6 ou shift F1 à F6 sur Casio, shift puis une touche 1 à 0, et les parenthèses sur Nspire et Numworks). Les mots clefs apparaissent en couleur et il y a une aide à la saisie des parenthèses et crochets. La touche curseur vers le bas permet de compléter une commande dont on a saisi le début ou/et d’afficher de l’aide et exemples sur cette commande.

Pour bloquer le clavier en mode alphabétique : sur les Casio, taper sur la touche F5, sur les Numworks taper 2 fois la touche alpha. Pour quitter le mode alphabétique, taper sur la touche alpha.

Exemples :

8.3  Programmation

Pour saisir un programme, depuis le shell KhiCAS, taper sur EXIT (Casio) ou esc (Nspire) ou Back (Numworks). Vous pouvez choisir entre syntaxe Xcas ou Python (shift SETUP sur les Casio). Le parser détecte aussi automatiquement la syntaxe, par exemple si vous commencez à taper def en début de programme. Les menus rapides F1 à F3 sur Casio et shift 1 à shift 3 sur TI/Numworks contiennent les structures de programmation les plus utiles.

Pour tester la syntaxe d’un programme, taper sur la touche de validation (EXE, enter ou OK). Utiliser shift-EXE pour passer à la ligne sur Casio, ret sur TI Nspire, ou EXE sur Numworks. S’il n’y a pas d’erreurs, pensez à sauvegarder votre programme avant de l’exécuter pour ne pas risquer de tout perdre en cas de crash.

Puis vous pouvez taper EXIT, esc ou Back pour revenir au shell, d’où vous pouvez exécuter le programme. On peut interrompre l’exécution d’un programme en appuyant sur la touche AC/ON ou ON ou Back selon la calculatrice. Si un programme ne fonctionne pas comme prévu, la commande debug() permet de l’exécuter en pas à pas en affichant les variables locales et permet de le mettre au point.

9  Quelques références

Index

A  La moyenne arithmético-géométrique.

A.1  Définition et convergence

Soient aa et bb deux réels positifs, on définit les 2 suites u 0=a,v 0=b,u n+1=u n+v n2,v n+1=u nv n(15) u_0=a, v_0=b, \quad u_{n+1}=\frac{u_n+v_n}{2}, v_{n+1}=\sqrt{u_nv_n} \qquad (15) On va montrer que ces 2 suites sont adjacentes et convergent donc vers une limite commune notée M(a,b)M(a,b) et il se trouve que la convergence est très rapide, en raison de l’identité : u n+1v n+1=12(u nv n) 2=12(u n+v n) 2(u nv n) 2(16) u_{n+1}-v_{n+1}=\frac{1}{2}(\sqrt{u_n}-\sqrt{v_n})^2 =\frac{1}{2(\sqrt{u_n}+\sqrt{v_n})^2}(u_n-v_n)^2 \qquad (16) la convergence est quadratique.

On suppose dans la suite que aba\geq b sans changer la généralité puisque échanger aa et bb ne change pas la valeur de u nu_n et v nv_n pour n>0n&gt;0. On a alors u nv nu_n \geq v_n (d’après (16) pour n>0n&gt;0) et u n+1u nu_{n+1} \leq u_n car u n+1u n=12(v nu n)0 u_{n+1}-u_n=\frac{1}{2}(v_n-u_{n}) \leq 0 et v n+1=u nv nv nv n=v nv_{n+1}=\sqrt{u_nv_n} \geq \sqrt{v_nv_n}=v_n. Donc (u n)(u_n) est décroissante minorée (par v 0v_0), (v n)(v_n) est croissante majorée (par u 0u_0), ces 2 suites sont convergentes et comme u n+1=u n+v n2u_{n+1}=\frac{u_n+v_n}{2}, elles convergent vers la même limite ll qui dépend de aa et bb et que l’on note M(a,b)M(a,b). On remarque aussi que M(a,b)=bM(a/b,1)=aM(1,b/a)M(a,b)=bM(a/b,1)=aM(1,b/a).

Précisons maintenant la vitesse de convergence lorsque ab>0a \geq b&gt;0. On va commencer par estimer le nombre d’itérations nécessaires pour que u nu_n et v nv_n soient du même ordre de grandeur. Pour cela, on utilise la majoration ln(u n+1)ln(v n+1)ln(u n)ln(v n+1)=12(ln(u n)ln(v n)) \ln(u_{n+1})-\ln(v_{n+1}) \leq \ln(u_{n})-\ln(v_{n+1}) = \frac{1}{2} (\ln(u_{n})-\ln(v_{n})) donc lnu nv n=ln(u n)ln(v n)12 n(ln(a)ln(b))=12 nlnab \ln \frac{u_n}{v_n} = \ln(u_n)-\ln(v_n) \leq \frac{1}{2^n} (\ln(a)-\ln(b)) = \frac{1}{2^n} \ln \frac{a}{b} Donc si nln(ln(a/b)/m)ln(2)n \geq \frac{\ln( \ln(a/b)/m)}{\ln(2)} alors lnu nv nm\ln \frac{u_n}{v_n} \leq m (par exemple, on peut prendre m=0.1m=0.1 pour avoir u n/v n[1,e 0.1])u_n/v_n \in [1,e^{0.1}]). Le nombre minimum d’itérations n 0n_0 est proportionnel au log du log du rapport a/ba/b. Ensuite on est ramené à étudier la convergence de la suite arithmético-géométrique de premiers termes a=u n 0a=u_{n_0} et b=v n 0b=v_{n_0} et même en tenant compte de M(a,b)=aM(1,b/a)M(a,b)=aM(1,b/a) à a=1a=1 et b=v n/u nb=v_n/u_n donc 0ab1e 0.10\leq a-b \leq 1-e^{-0.1}. Alors l’équation (16) entraine u n+1v n+118(u nv n) 2 u_{n+1}-v_{n+1} \leq \frac{1}{8}(u_n-v_n)^2 puis (par récurrence) 0u nv n18 2 n1(ab) 2 n 0 \leq u_n-v_n \leq \frac{1}{8^{2^n-1}}(a-b)^{2^n} Donc comme M(a,b)M(a,b) est compris entre v nv_n et u nu_n, l’erreur relative sur la limite commune est inférieure à une précision donnée ε\epsilon au bout d’un nombre d’itérations proportionnel au ln(ln(1/ε))\ln(\ln(1/\epsilon)).

Typiquement dans la suite, on souhaitera calculer M(1,b)M(1,b) avec bb de l’ordre de 2 n2^{-n} en déterminant nn chiffres significatifs, il faudra alors O(ln(n))O(\ln(n)) itérations pour se ramener à M(1,b)M(1,b) avec b[e 0.1,1]b\in [e^{-0.1},1] puis O(ln(n))O(\ln(n)) itérations pour avoir la limite avec nn chiffres significatifs.

Le cas complexe
On suppose maintenant que a,ba, b \in \mathbb{C} avec (a)>0,(b)>0\Re(a)&gt;0, \Re(b)&gt;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 (u n)>0\Re(u_n)&gt;0 ; de plus (v n)>0\Re(v_n) &gt; 0 car par définition de la racine carrée (v n)0\Re(v_n)\geq 0 et est de plus non nul car le produit de deux complexes d’arguments dans ]π/2,π/2[]-\pi/2,\pi/2[ ne peut pas être un réel négatif. On en déduit que arg(u n+1)=arg(u n+v n)\arg(u_{n+1})=\arg(u_n+v_n) se trouve dans l’intervalle de bornes arg(u n)\arg(u_n) et arg(v n)\arg(v_n) et que arg(v n+1)=12(arg(u n)+arg(v n))\arg(v_{n+1})=\frac{1}{2}(\arg(u_n)+\arg(v_n)) donc |arg(u n+1arg(v n+1)|12|arg(u n)arg(v n)| | \arg(u_{n+1}-\arg(v_{n+1}) | \leq \frac{1}{2}|\arg(u_n)-\arg(v_n)| Après nn itérations, on a |arg(u n)arg(v n)|π2 n |\arg(u_n)-\arg(v_n)| \leq \frac{\pi}{2^n} Après quelques itérations, u nu_n et v nv_n seront donc presque alignés. Faisons 4 itérations. On peut factoriser par exemple v nv_n et on est ramené à l’étude de la suite de termes initiaux a=u n/v na=u_n/v_n d’argument arg(u n)arg(v n)\arg(u_n)-\arg(v_n) petit (inférieur en valeur absolue à π/16\pi/16) et b=1b=1. On suppose donc dans la suite que |arg(u nv n)|π/162 n |\arg(\frac{u_n}{v_n})| \leq \frac{\pi/16}{2^n} Étude du module
On a : u n+1v n+1=12(u nv n+1u nv n) \frac{u_{n+1}}{v_{n+1}}= \frac{1}{2}\left(\sqrt{\frac{u_{n}}{v_{n}}}+\frac{1}{\sqrt{\frac{u_{n}}{v_{n}}}}\right) Posons u nv n=ρ ne iθ n \frac{u_{n}}{v_{n}}=\rho_n e^{i\theta_n}, on a : |u n+1v n+1| = 12|ρ ne iθ n/2+1ρ ne iθ n/2| = 12|(ρ n+1ρ n)cosθ n2+i(ρ n1ρ n)sinθ n2| = 12(ρ n+1ρ n) 2cos 2θ n2+(ρ n1ρ n) 2sin 2θ n2 = 12ρ n+1ρ n+2cosθ n \begin{matrix} |\frac{u_{n+1}}{v_{n+1}}| &= & \frac{1}{2}\left|\sqrt{\rho_n} e^{i\theta_n/2} + \frac{1}{\sqrt{\rho_n}} e^{-i\theta_n/2} \right| \\ &=& \frac{1}{2} \left| (\sqrt{\rho_n}+ \frac{1}{\sqrt{\rho_n}})\cos\frac{\theta_n}{2} + i (\sqrt{\rho_n}- \frac{1}{\sqrt{\rho_n}})\sin\frac{\theta_n}{2} \right| \\ &=& \frac{1}{2} \sqrt{ (\sqrt{\rho_n}+ \frac{1}{\sqrt{\rho_n}})^2\cos^2\frac{\theta_n}{2} + (\sqrt{\rho_n}- \frac{1}{\sqrt{\rho_n}})^2\sin^2\frac{\theta_n}{2} } \\ &=& \frac{1}{2} \sqrt{ \rho_n+ \frac{1}{\rho_n} +2\cos \theta_n } \end{matrix} Si ρ\rho désigne le max de ρ n\rho_n et 1/ρ n1/\rho_n, on a alors la majoration |u n+1v n+1|12ρ+ρ+2ρ=ρ |\frac{u_{n+1}}{v_{n+1}}| \leq \frac{1}{2} \sqrt{ \rho + \rho + 2 \rho } = \sqrt{\rho} donc en prenant les logarithmes lnρ n+112lnρ=12|lnρ n|(17) \ln \rho_{n+1} \leq \frac{1}{2} \ln \rho=\frac{1}{2} |\ln \rho_n| \qquad (17) On rappelle qu’on a la majoration |arg(u nv n)|=|θ n|π/162 n12 n+1 |\arg(\frac{u_n}{v_n})| = |\theta_n| \leq \frac{\pi/16}{2^n} \leq \frac{1}{2^{n+1}} qui va nous donner la minoration de ρ n+1\rho_{n+1} ρ n+1=|u n+1v n+1| = 12ρ n+1ρ n+22(1cosθ n) = 12ρ n+1ρ n+24sin 2(θ n2) 12ρ n+1ρ n+2θ n 2 12ρ n+1ρ n+2×1θ n 2ρ n+1ρ n+2 121ρ+1ρ+21ρ×1θ n 24 1ρ1θ n 24 1ρ114×2 2n+2 \begin{matrix} \rho_{n+1}=|\frac{u_{n+1}}{v_{n+1}}| & = & \frac{1}{2} \sqrt{ \rho_n+ \frac{1}{\rho_n} +2 - 2 (1-\cos \theta_n) } \\ & = & \frac{1}{2} \sqrt{ \rho_n+ \frac{1}{\rho_n} +2 - 4 \sin^2 (\frac{\theta_n}{2}) } \\ & \geq & \frac{1}{2} \sqrt{ \rho_n+ \frac{1}{\rho_n} +2 - \theta_n^2} \\ & \geq & \frac{1}{2} \sqrt{ \rho_n+ \frac{1}{\rho_n} +2} \times \sqrt{1 - \frac{\theta_n^2}{\rho_n+ \frac{1}{\rho_n} +2}} \\ & \geq & \frac{1}{2} \sqrt{ \frac{1}{\rho} + \frac{1}{\rho} +2\frac{1}{\rho}} \times \sqrt{1 - \frac{\theta_n^2}{4}} \\ & \geq & \frac{1}{\sqrt{\rho}} \sqrt{1 - \frac{\theta_n^2}{4}} \\ & \geq & \frac{1}{\sqrt{\rho}} \sqrt{1 - \frac{1} {4 \times 2^{2n+2}}} \end{matrix} en prenant les log et en minorant ln(1x)\ln(1-x) par 2x-2x lnρ n+112(|lnρ n|+ln(114×2 2n+2))12(|lnρ n|+12 2n+3) \ln \rho_{n+1} \geq \frac{1}{2} (-|\ln \rho_n|+\ln(1 -\frac{1} {4 \times 2^{2n+2}} )) \geq -\frac{1}{2} (|\ln \rho_n|+\frac{1} {2^{2n+3}} ) Finalement avec (17) |lnρ n+1|12(|lnρ n|+12 2n+3) |\ln \rho_{n+1}| \leq \frac{1}{2} (|\ln \rho_n|+\frac{1}{2^{2n+3}} ) On en déduit |lnρ n|12 nlnρ 0+12 n+3+...+12 2n+1+12 2n+2=12 nlnρ 0+12 n+2 |\ln \rho_n| \leq \frac{1}{2^n} \ln \rho_0 + \frac{1}{2^{n+3}} + ... + \frac{1}{2^{2n+1}} + \frac{1}{2^{2n+2}} = \frac{1}{2^n} \ln \rho_0 + \frac{1}{2^{n+2}} La convergence du ln(u n/v n)\ln(u_n/v_n) vers 0 est donc géométrique, donc u nu_n et v nv_n convergent quadratiquement.

A.2  Lien avec les intégrales elliptiques

Le calcul de la limite commune des suites u nu_n et v nv_n en fonction de aa et bb 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 aa et bb et qui est invariante par la transformation u n,v nu n+1,v n+1u_n,v_n \rightarrow u_{n+1},v_{n+1} (15) I(a,b)= +dt(a 2+t 2)(b 2+t 2) I(a,b)=\int_{-\infty}^{+\infty} \frac{dt} {\sqrt{(a^2+t^2)(b^2+t^2)}} On a en effet I(a+b2,ab)= +du((a+b2) 2+u 2)(ab+u 2) I(\frac{a+b}{2},\sqrt{ab}) = \int_{-\infty}^{+\infty} \frac{du}{\sqrt{((\frac{a+b}{2})^2+u^2)(ab+u^2)}} On pose alors u=12(tabt),t>0 u=\frac{1}{2} (t-\frac{ab}{t}), \quad t&gt;0 tut \rightarrow u est une bijection croissante de t]0,+[t\in]0,+\infty[ vers u],+[u \in ]-\infty,+\infty[, donc I(a+b2,ab) = 0 +dt/2(1+ab/t 2)((a+b2) 2+1/4(tab/t) 2)(ab+1/4(tab/t) 2) = 2 0 +dt(a 2+t 2)(b 2+t 2)=I(a,b) \begin{matrix} I(\frac{a+b}{2},\sqrt{ab}) &=& \int_{0}^{+\infty} \frac{dt/2(1+ab/t^2)}{\sqrt{((\frac{a+b}{2})^2+1/4(t-ab/t)^2)(ab+1/4(t-ab/t)^2)}}\\ &=& 2 \int_{0}^{+\infty} \frac{dt}{\sqrt{(a^2+t^2)(b^2+t^2)}} = I(a,b) \end{matrix} On note au passage que II est définie si a,ba,b \in \mathbb{C} vérifient (a)>0,(b)>0\Re(a)&gt;0, \Re(b)&gt;0, on peut montrer que la relation ci-dessus s’étend (par holomorphie).

Lorsque a=b=la=b=l (par exemple lorsqu’on est à la limite), le calcul de I(l,l)I(l,l) est explicite I(l,l)= +dt(l 2+t 2)=πl I(l,l)=\int_{-\infty}^{+\infty} \frac{dt}{(l^2+t^2)} = \frac{\pi}{l} donc I(a,b)=I(M(a,b),M(a,b))=πM(a,b) I(a,b)=I(M(a,b),M(a,b))=\frac{\pi}{M(a,b)} On peut transformer I(a,b)I(a,b) en posant t=but=bu I(a,b)=2 0 +du(a 2+b 2u 2)(1+u 2)=2a 0 +du(1+(b/a) 2u 2)(1+u 2) I(a,b)=2\int_{0}^{+\infty} \frac{du}{\sqrt{(a^2+b^2u^2)(1+u^2)}} = \frac{2}{a} \int_{0}^{+\infty} \frac{du}{\sqrt{(1+(b/a)^2u^2)(1+u^2)}} Puis en posant u=tan(x)u=\tan(x) (du=(1+u 2)dxdu=(1+u^2) dx) I(a,b)=2a 0 π21+tan(x) 21+(b/a) 2tan(x) 2dx I(a,b)=\frac{2}{a} \int_0^{\frac{\pi}{2}} \sqrt{\frac{1+\tan(x)^2}{1+(b/a)^2\tan(x)^2}} \ dx et enfin en posant tan 2(x)=sin(x) 21sin(x) 2\tan^2(x)=\frac{\sin(x)^2}{1-\sin(x)^2} I(a,b)=2a 0 π211(1b 2a 2)sin(x) 2dx I(a,b)= \frac{2}{a} \int_0^{\frac{\pi}{2}} \sqrt{ \frac{1}{1-(1-\frac{b^2}{a^2})\sin(x)^2} } \ dx Si on définit pour m<1m&lt;1 K(m)= 0 π2dx1msin(x) 2 K(m)=\int_0^{\frac{\pi}{2}} \frac{dx}{\sqrt{1-m \sin(x)^2}} alors on peut calculer KK en fonction de II, en posant m=1b 2/a 2m=1-b^2/a^2 soit b 2/a 2=1mb^2/a^2=1-m K(m)=a2I(a,a1m)=a2πM(a,a1m)=π2M(1,1m) K(m)=\frac{a}{2} I(a,a\sqrt{1-m})=\frac{a}{2}\frac{\pi}{M(a,a\sqrt{1-m})} =\frac{\pi}{2M(1,\sqrt{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)= 0 π2dx1msin(x) 2=π2M(1,1m)(18) K(m)=\int_0^{\frac{\pi}{2}} \frac{dx}{\sqrt{1-m \sin(x)^2}}= \frac{\pi}{2M(1,\sqrt{1-m})} \qquad (18) Dans l’autre sens, pour xx et yy positifs K((xyx+y) 2)=π2M(1,1(xyx+y) 2)=π2M(1,2x+yxy)=π22x+yM(x+y2,xy)=π4x+yM(x,y) K( (\frac{x-y}{x+y})^2 )= \frac{\pi}{2M(1,\sqrt{1-(\frac{x-y}{x+y})^2})} = \frac{\pi}{2M(1,\frac{2}{x+y}\sqrt{xy})} = \frac{\pi}{2 \frac{2}{x+y} M(\frac{x+y}{2},\sqrt{xy}) } = \frac{\pi}{4} \frac{x+y}{M(x,y)} et finalement M(x,y)=π4x+yK((xyx+y) 2) M(x,y)=\frac{\pi}{4} \frac{x+y}{ K\left( (\frac{x-y}{x+y}\right)^2 )}

A.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)K(m) lorsque mm tend vers 1. Plus précisément, on va poser 1m=k 21-m=k^2 avec k]0,1]k \in ]0,1], donc K(m)= 0 π2dx1(1k 2)sin(x) 2= 0 π2dy1(1k 2)cos(y) 2 K(m)= \int_0^{\frac{\pi}{2}} \frac{dx}{\sqrt{1-(1-k^2) \sin(x)^2}} =\int_0^{\frac{\pi}{2}} \frac{dy}{\sqrt{1-(1-k^2) \cos(y)^2}} en posant y=π/2xy=\pi/2-x, et K(m)= 0 π2dysin(y) 2+k 2cos(y) 2 K(m)=\int_0^{\frac{\pi}{2}} \frac{dy}{\sqrt{\sin(y)^2+k^2 \cos(y)^2}} la singularité de l’intégrale pour kk proche de 0 apparait lorsque yy est proche de 0. Si on effectue un développement de Taylor en y=0y=0, on trouve sin(y) 2+k 2cos(y) 2=k 2+(1k 2)y 2+O(y 4) \sin(y)^2+k^2 \cos(y)^2 = k^2 + (1-k^2) y^2 + O(y^4) Il est donc naturel de comparer K(m)K(m) à l’intégrale J= 0 π2dyk 2+(1k 2)y 2 J=\int_0^{\frac{\pi}{2}} \frac{dy}{\sqrt{k^2 + (1-k^2) y^2}} qui se calcule en faisant par exemple le changement de variables y=k1k 2sinh(t) y=\frac{k}{\sqrt{1-k^2}} \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=11k 2(ln(πk)+ln(12(1k 2+4k 2π 2+1k 2)))(19) J= \frac{1}{\sqrt{1-k^{2}}} \left( \ln\left(\frac{\pi}{k}\right) + \ln\left( \frac{1}{2} \left(\sqrt{ 1-k^{2} +4 \frac{k^{2}}{\pi^2}} +\sqrt{1-k^{2}} \right) \right) \right) \qquad (19) et on peut calculer le développement asymptotique de JJ en 0

series(J,k=0,5,1)

qui renvoie : J=ln(πk)+O((1ln(k)) 5) J =\ln\left(\frac{\pi}{k}\right) +O( \left(\frac{-1}{\ln(k)}\right)^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 O˜\tilde{O} peut contenir des logarithmes) (1π 2+ln(π)ln(k)12)k 2+O˜(k 4) \left(\frac{1}{\pi^2} + \frac{\ln(\pi)-\ln(k)-1}{2}\right) k^{2} + \tilde{O}(k^4) donc J=ln(k)+ln(π)+(1π 2+ln(π)ln(k)12)k 2+O˜(k 4)(20) J=-\ln(k)+\ln(\pi)+\left(\frac{1}{\pi^2} + \frac{\ln(\pi)-\ln(k)-1}{2}\right) k^{2} + \tilde{O}(k^4) \qquad (20) Examinons maintenant KJK-J, il n’a plus de singularité en y=0y=0, et il admet une limite lorsque k0k \rightarrow 0, obtenue en remplacant kk par 0 (KJ) |k=0= 0 π2(1sin(y)1y)dy=[ln(tan(y2))ln(y)] 0 π2=ln(4π) (K-J)_{|k=0} = \int_0^{\frac{\pi}{2}} \left(\frac{1}{\sin(y)}-\frac{1}{y}\right) \ dy = \left[\ln\left(\tan\left(\frac{y}{2}\right)\right) - \ln(y) \right]_0^{\frac{\pi}{2}} = \ln(\frac{4}{\pi}) D’où pour KK K k0=ln(4k)+O((1ln(k)) 5) K_{k \rightarrow 0} = \ln\left(\frac{4}{k}\right) + O( \left(\frac{-1}{\ln(k)}\right)^5) Pour préciser la partie du développement de KK en puissances de kk, nous allons majorer KJln(4/π)K-J-\ln(4/\pi), puis Jln(π/k)J-\ln(\pi/k). Posons A=sin(y) 2+k 2cos(y) 2,B=y 2+(1y 2)k 2 A=\sin(y)^2+k^2 \cos(y)^2, \quad B=y^2+(1-y^2)k^2 Majoration de KJln(4/π)K-J-\ln(4/\pi)
L’intégrand de la différence KJln(4π)K-J-\ln(\frac{4}{\pi}) est 1A1B(1sin(y)1y) = BAABysin(y)ysin(y) = BAAB(A+B)ysin(y)ysin(y) = (y 2sin(y) 2)(1k 2)AB(A+B)ysin(y)ysin(y)(21) \begin{matrix} \frac{1}{\sqrt{A}} - \frac{1}{\sqrt{B}} - \left( \frac{1}{\sin(y)}-\frac{1}{y} \right) &= & \frac{\sqrt{B}-\sqrt{A}}{\sqrt{A} \sqrt{B}} - \frac{y-\sin(y)}{y\sin(y)} \\ &= & \frac{B-A}{\sqrt{A} \sqrt{B} (\sqrt{A}+\sqrt{B})} - \frac{y-\sin(y)}{y\sin(y)} \\ &=& \frac{(y^2-\sin(y)^2)(1-k^2)}{\sqrt{A} \sqrt{B} (\sqrt{A}+\sqrt{B})} - \frac{y-\sin(y)}{y\sin(y)} \end{matrix} \qquad (21) Soit KJln(4π)= 0 π2(ysin(y))[(1k 2)ysin(y)(y+sin(y))AB(A+B)]AB(A+B)ysin(y)(22) K-J-\ln(\frac{4}{\pi})= \int_0^{\frac{\pi}{2}} \frac{(y-\sin(y))[(1-k^2)y \sin(y)(y+\sin(y))-\sqrt{AB}(\sqrt{A}+\sqrt{B})]} {\sqrt{A} \sqrt{B} (\sqrt{A}+\sqrt{B})y\sin(y)} \qquad (22) On décompose l’intégrale en 2 parties [0,k][0,k] et [k,π/2][k,\pi/2]. Sur [0,k][0,k] on utilise (21), on majore chaque terme séparément et on minore AA et BB par A=k 2+(1k 2)sin(y) 2k 2,B=k 2+(1k 2)y 2k 2 A=k^2+(1-k^2)\sin(y)^2 \geq k^2, \quad B=k^2+(1-k^2)y^2 \geq k^2 Donc | 0 k| 0 k|BA|2k 3dy+ 0 k(1sin(y)1y)dy 0 ky 2sin(y) 22k 3dy+ln(tan(k2))ln(k2) 13k 3+12k+14sin(2k)2k 3+ln(sin(k2))ln(k2)ln(cos(k2)) 13k 3+12k+14(2k8k 36+32k 55!2k 3ln(cos(k2)) k 230ln(112!(k2) 2) k 230+k 24 \begin{matrix} | \int_0^{k} | &\leq &\int_0^k \frac{|B-A|}{2k^3} \ dy + \int_0^k ( \frac{1}{\sin(y)}-\frac{1}{y} ) \ dy \\ &\leq& \int_0^k \frac{y^2-\sin(y)^2}{2k^3} \ dy + \ln (\tan(\frac{k}{2})) -\ln(\frac{k}{2}) \\ &\leq & \frac{\frac{1}{3} k^{3}+\frac{-1}{2} k+\frac{1}{4} \sin(2 k)}{2 k^{3}} + \ln (\sin(\frac{k}{2})) -\ln(\frac{k}{2}) - \ln (\cos(\frac{k}{2})) \\ &\leq & \frac{\frac{1}{3} k^{3}+\frac{-1}{2} k+\frac{1}{4} (2k-\frac{8k^3}{6}+\frac{32k^5}{5!}}{2 k^{3}} - \ln (\cos(\frac{k}{2})) \\ &\leq & \frac{k^2}{30}- \ln (1- \frac{1}{2!}\left(\frac{k}{2}\right)^2) \\ &\leq & \frac{k^2}{30} +\frac{k^2}{4} \end{matrix} Sur [k,π/2][k,\pi/2], on utilise (22) et on minore AA et BB par A=sin(y) 2+k 2cos(y) 2sin(y) 2,B=y 2+(1y 2)k 2y 2 A=\sin(y)^2+k^2 \cos(y)^2 \geq \sin(y)^2, \quad B=y^2+(1-y^2)k^2 \geq y^2 on obtient | k π2| k π2(ysin(y))|C|ysin(y)(y+sin(y)), | \int_k^{\frac{\pi}{2}} | \leq \int_k^{\frac{\pi}{2}} \frac{(y-\sin(y))|C|} {y \sin(y) (y+\sin(y))} , où : C = (1k 2)ysin(y)(y+sin(y))AB+BA = A(By)B(Asin(y))AyBsin(y)+(1k 2)ysin(y)(y+sin(y)) = A(By)B(Asin(y))k 2(y+sin(y)) \begin{matrix} C&=&(1-k^2)y \sin(y)(y+\sin(y))-A\sqrt{B}+B\sqrt{A} \\ &=& -A(\sqrt{B}-y)-B(\sqrt{A}-\sin(y)) -Ay-B\sin(y) + (1-k^2)y \sin(y)(y+\sin(y)) \\ &=& -A(\sqrt{B}-y)-B(\sqrt{A}-\sin(y)) - k^2(y+\sin(y)) \end{matrix} Donc |C| A(By)+B(Asin(y))+k 2(y+sin(y)) ABy 2B+y+BAsin(y) 2A+sin(y)+k 2(y+sin(y)) Ak 22y+Bk 22sin(y)+k 2(y+sin(y)) \begin{matrix} |C| &\leq& A(\sqrt{B}-y)+B(\sqrt{A}-\sin(y)) + k^2(y+\sin(y)) \\ &\leq& A \frac{B-y^2}{\sqrt{B}+y} + B \frac{A-\sin(y)^2}{\sqrt{A}+\sin(y)} + k^2(y+\sin(y)) \\ &\leq & A \frac{k^2}{2y} + B \frac{k^2}{2\sin(y)} + k^2(y+\sin(y)) \end{matrix} et | k π2| k π2(ysin(y))k 2(A2y+B2sin(y)+(y+sin(y)))ysin(y)(y+sin(y)) | \int_k^{\frac{\pi}{2}} | \leq \int_k^{\frac{\pi}{2} } \frac{(y-\sin(y))k^2(\frac{A}{2y} + \frac{B}{2\sin(y)} + (y+\sin(y))) } {y \sin(y) (y+\sin(y))} On peut majorer ysin(y)y 3/6y-\sin(y) \leq y^3/6, donc | k π2|k 26 k π2Ay2sin(y)(sin(y)+y)+By 2sin(y) 2(sin(y)+y)+y 2sin(y) | \int_k^{\frac{\pi}{2}} | \leq \frac{k^2}{6} \int_k^{\frac{\pi}{2}} \frac{Ay}{2\sin(y) (\sin(y)+y)} + \frac{By^2}{\sin(y)^2(\sin(y)+y)} + \frac{y^2}{\sin(y)} On majore enfin AA et BB par 1, | k π2|k 26 k π2y2sin(y) 2+y 2sin(y) | \int_k^{\frac{\pi}{2}} | \leq \frac{k^2}{6} \int_k^{\frac{\pi}{2}} \frac{y}{2\sin(y)^2} + \frac{y^2}{\sin(y)} Le premier morceau se calcule par intégration par parties k 26 k π2y2sin(y) 2 = k 26([ytan(y)] k π/2+ k π21tan(y)) = k 26(ktan(k)+[ln(sin(y))] k π2) = k 26(ktan(k)ln(sin(k))) k 26(1ln(k)) \begin{matrix} \frac{k^2}{6} \int_k^{\frac{\pi}{2}} \frac{y}{2\sin(y)^2} &=& \frac{k^2}{6} \left( [-\frac{y}{\tan(y)}]_k^{\pi/2} + \int_k^{\frac{\pi}{2}} \frac{1}{\tan(y)} \right) \\ &=& \frac{k^2}{6} \left(\frac{k}{\tan(k)}+ [\ln(\sin(y))]_k^{\frac{\pi}{2}} \right)\\ &=& \frac{k^2}{6} \left(\frac{k}{\tan(k)}-\ln(\sin(k)) \right)\\ &\leq & \frac{k^2}{6}(1-\ln(k)) \end{matrix} Le deuxième morceau se majore en minorant sin(y)(2y)/π\sin(y)\geq (2y)/\pi k 26 k π2y 2sin(y)k 26 0 π2π2y=k 2π 396 \frac{k^2}{6} \int_k^{\frac{\pi}{2}} \frac{y^2}{\sin(y)} \leq \frac{k^2}{6} \int_0^{\frac{\pi}{2}} \frac{\pi}{2} y = \frac{k^2\pi^3}{96} Finalement |KJln(4π)|k 2(16ln(k)+π 396+16+130+14) |K-J-\ln(\frac{4}{\pi})| \leq k^2 \left( -\frac{1}{6} \ln(k) + \frac{\pi^3}{96} + \frac{1}{6} + \frac{1}{30}+ \frac{1}{4} \right) JJ est donné en (19).

Majoration de Jln(π/k)J-ln(\pi/k)
On a |Jln(πk)|=|(11k 21)ln(πk)+11k 2ln(12(1k 2+4k 2π 2+1k 2))| |J - \ln\left(\frac{\pi}{k}\right)| = \left| (\frac{1}{\sqrt{1-k^2}}-1) \ln\left(\frac{\pi}{k}\right) + \frac{1}{\sqrt{1-k^2}} \ln\left( \frac{1}{2} \left(\sqrt{ 1-k^{2} +4 \frac{k^{2}}{\pi^2}} +\sqrt{1-k^{2}} \right) \right) \right| et on va majorer la valeur absolue de chaque terme de la somme. Pour k1/2k\leq 1/2, on a 11k 21=k 21k 2+1k 2k 23/4+3/2 \frac{1}{\sqrt{1-k^2}}-1=\frac{k^2}{\sqrt{1-k^2}+1-k^2} \leq \frac{k^2}{3/4+\sqrt{3}/2} Pour le second terme, on majore le facteur 11k 2\frac{1}{\sqrt{1-k^2}} par 23\frac{2}{\sqrt{3}}, l’argument du logarithme est inférieur à 1 et supérieur à 12(1k 22+1k 2(14π 2)2)=1k 2(11π 2)>1k 2 \frac{1}{2}(1 - \frac{k^2}{2} +1- \frac{k^2(1-\frac{4}{\pi^2})}{2}) = 1 - k^2 ( 1-\frac{1}{\pi^2}) &gt; 1-k^2 donc le logarithme en valeur absolue est inférieur à 2k 2 2 k^2 donc, pour k1/2k\leq 1/2, |Jln(πk)|k 23/4+3/2ln(πk)+k 243 |J-\ln\left(\frac{\pi}{k}\right)| \leq \frac{k^2}{3/4+\sqrt{3}/2} \ln\left(\frac{\pi}{k}\right) + k^2 \frac{4}{\sqrt{3}} Finalement, pour k<1/2k&lt;1/2 |Kln(4k)|k 2(lnπ3/4+3/2+43+π 396+920(13/4+3/2+16)ln(k))(23) |K-\ln\left(\frac{4}{k}\right) | \leq k^2 \left( \frac{\ln \pi}{3/4+\sqrt{3}/2} + \frac{4}{\sqrt{3} } + \frac{\pi^3}{96} + \frac{9}{20} - (\frac{1}{3/4+\sqrt{3}/2}+\frac{1}{6}) \ln(k) \right) \qquad (23) que l’on peut réécrire |π2M(1,k)ln(4k)|k 2(3.80.8ln(k))(24) |\frac{\pi}{2M(1,k)}-\ln\left(\frac{4}{k}\right) | \leq k^2(3.8-0.8\ln(k)) \qquad (24) La formule (24) permet de calculer le logarithme d’un réel positif avec (presque) nn bits lorsque k2 n/2k \leq 2^{-n/2} (ce à quoi on peut toujours se ramener en calculant le logarithme d’une puissance 2 m2^m-ième de xx ou le logarithme de 2 mx2^{m}x, en calculant au préalable ln(2)\ln(2)). Par exemple, prenons k=2 27k=2^{-27}, on trouve (en 8 itérations) M(1,2 27)=M 1=0.0781441403763M(1,2^-{27})=M_1=0.0781441403763. On a, avec une erreur inférieure à 19×2 54=1.1×10 1519 \times 2^{-54}=1.1\times 10^{-15} M(1,2 27)=M 1=π2ln(2 29)=π58ln(2), M(1,2^-{27})=M_1=\frac{\pi}{2\ln(2^{29})}=\frac{\pi}{58\ln(2)}, On peut donc déduire une valeur approchée de π\pi si on connait la valeur approchée de ln(2)\ln(2) et réciproquement. Si on veut calculer les deux simultanément, comme les relations entre ln\ln et π\pi seront des équations homogènes, on est obligé d’introduire une autre relation. Par exemple pour calculer une valeur approchée de π\pi on calcule la différence ln(2 29+1)ln(2 29)\ln(2^{29}+1)-\ln(2^{29}) dont on connait le développement au premier ordre, et on applique la formule de la moyenne arithmético-géométrique. Il faut faire attention à la perte de précision lorsqu’on fait la différence des deux logarithmes qui sont très proches, ainsi on va perdre une trentaine de bits, il faut grosso modo calculer les moyennes arithmético-géométrique avec 2 fois plus de chiffres significatifs.

L’intérêt de cet algorithme apparait lorsqu’on veut calculer le logarithme avec beaucoup de précision, en raison de la convergence quadratique de la moyenne arithmético-géométrique (qui est nettement meilleure que la convergence linéaire pour les développements en série, ou logarithmiquement meilleure pour l’exponentielle), par contre elle n’est pas performante si on ne veut qu’une dizaine de chiffres significatifs. On peut alors calculer les autres fonctions transcendantes usuelles, telle l’exponentielle, à partir du logarithme, ou les fonctions trigonométriques inverses (en utilisant des complexes) et directes.

On trouvera dans Brent-Zimmermann quelques considérations permettant d’améliorer les constantes dans les temps de calcul par rapport à cette méthode (cela nécessite d’introduire des fonctions spéciales θ\theta) et d’autres formules pour calculer π\pi.

  

Retour à la page principale de Giac/Xcas.