\documentclass{article} \usepackage[latin1]{inputenc} \usepackage[T1]{fontenc} \usepackage{amsmath,amsfonts,amssymb} \usepackage[francais]{babel} \usepackage{times} \usepackage{xspace} \usepackage{makeidx} \usepackage{ifpdf} \ifpdf \usepackage[colorlinks,pdftex]{hyperref} \else \usepackage[ps2pdf,breaklinks=true,colorlinks=true,linkcolor=red,citecolor=green]{hyperref} \fi \newcommand{\bs}{\symbol{92}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\Q}{\mathbb{Q}} \newcommand{\N}{\mathbb{N}} \newcommand{\R}{\mathbb{R}} \newcommand{\C}{\mathbb{C}} \newtheorem{theorem}{Théorème} \newtheorem{prop}{Proposition}[theorem] \newtheorem{defn}{Définition} \makeindex %HEVEA\renewcommand{\cuttingunit}{subsection} \begin{document} \title{Mat 249} \author{Bernard.Parisse@ujf-grenoble.fr} \date{2009/10} \maketitle \tableofcontents \section{Présentation du module} %\verb|^|\\ Dans ce module, on introduira moins de notions nouvelles que dans d'autres modules de math\'ematiques, par contre on insistera sur le calcul effectif, si possible efficace, et sur le controle de la pr\'ecision des r\'esultats, ceci explique la part importante consacr\'ee aux TP (18h de cours en 12 séances, 18h de TD en 12 séances et 24h de TP). On pr\'esentera par exemple des m\'ethodes de calcul des fonctions usuelles (racine carr\'ee, trigonom\'etriques, ...), il s'agira non seulement de savoir calculer une valeur num\'erique, mais aussi de pouvoir majorer l'\'ecart entre la valeur trouv\'ee et la valeur exacte, en utilisant des th\'eor\`emes du cours. Les calculs se feront dans la mesure du possible sur ordinateur ou sur calculatrices. Les th\`emes abord\'es seront~: \begin{enumerate} \item calcul exact et approch\'e, repr\'esentation des donn\'ees \item Suites r\'ecurrentes, m\'ethode du point fixe, de Newton \item Interpolation polyn\^omiale (\'evaluation, interpolation de Lagrange) \item S\'eries de Taylor et approximation des fonctions usuelles \item Arithm\'etique des polyn\^omes (PGCD, B\'ezout, factorisation, d\'ecomposition en \'el\'ements simples) \item Int\'egration \item Alg\`ebre lin\'eaire \end{enumerate} L'évaluation se fait sur~: \begin{itemize} \item 1/4~: un DS à mi-semestre \item 1/4~: certains compte-rendus de TP (à rédiger seul ou en binome), \item 1/2~: l'examen final \end{itemize} Les calculatrices et les netbooks de taille d'écran plus petits que 10 pouces sont autorisées au DS et à l'examen final (pr\^et possible de netbooks pour le semestre). \section{Représentation des nombres et autres donn\'ees, calcul exact/approch\'e} Résumé:\\ Types de base~: entier machine, entier long, flottant machine et multipr\'ecision (Base 2, base 10, BCD). \\ Types compos\'es~: 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 $\N$ et relatifs $\Z$, les rationnels $\Q$, les réels $\R$ et les complexes $\C$. Sur ordinateur, on peut représenter ces nombres de manière exacte dans certains cas, approchée dans d'autres. \subsection{Entiers courts et longs} \begin{prop} {\bf Division euclidienne de deux entiers~:} si $a$ et $b$ sont deux entiers, $a \geq 0, b>0$, il existe un unique couple $(q,r)$ tel que \index{division euclidienne} \[ a = bq +r , \quad r \in [0, |b|[ \] \end{prop} Preuve~: On prend pour $q$ le plus grand entier tel que $a-bq>0$. La division euclidienne permet d'écrire un nombre entier, en utilisant une base $b$ et des caractères pour représenter les entiers entre 0 et $b-1$. Nous \'ecrivons les nombres entiers en base $b=10$ avec comme caractères les chiffres de 0 à 9. Les ordinateurs utilisent des circuits binaires pour stocker les informations, il est donc naturel d'y travailler en base 2 en utilisant comme caractères 0 et 1 ou en base 16 en utilisant comme caractères les chiffres de 0 à 9 et les lettres de A à F. En g\'en\'eral, pour trouver l'\'ecriture d'un nombre en base $b$ (par exemple $b=2$), on effectue des divisions euclidienne successives par $b$ du nombre puis de ses quotients successifs jusqu'\`a ce que le quotient soit 0 et on accolle les restes obtenus (premier reste \`a droite, dernier reste \`a gauche). Inversement, pour retrouver un entier $d$ \`a partir de son \'ecriture $d_n...d_0$, on traduit les divisions euclidiennes successives en \begin{eqnarray*} 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{eqnarray*} Par exemple, vingt-cinq s'\'ecrit en base 16 \verb|19| car 25 divis\'e par 16 donne quotient 1, reste 9. En base 2, on trouverait \verb|00011001| car $25=2^4+2^3+1$. On peut effectuer les op\'erations arithm\'etiques 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~: \begin{verbatim} 01001111 + 01101011 ---------- 10111010 \end{verbatim} {\bf 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\'eclin\'es en plusieurs variantes selon la taille et la possibilit\'e d'avoir un signe). Noter que la division de deux entiers $a$ et $b$ n'a pas la même signification que la division de deux réels, comme elle ne tomberait pas forcément juste, on calcule le quotient et le reste de la division euclidienne. Ces entiers machines permettent de représenter de manière exacte des petits entiers relatifs par exemple un entier machine sign\'e sur 4 octets est compris entre $[-2^{31},2^{31}-1]$. Selon le microprocesseur les 4 octets repr\'esentant l'entier sont stock\'es par adresse mémoire décroissante ou croissante (big ou little endian). Sur certains systèmes (dits BCD), on écrit les entiers en base 10, chaque chiffre occupant 4 bits (qui normalement sert à stocker un chiffre en base 16). Les microprocesseurs correspondants ont un flag leur permettant d'effectuer les opérations sur des nombres vu en représentation BCD (base 10) ou hexadécimale (base 16). Ces entiers machines permettent de faire très rapidement du calcul exact sur les entiers, mais \`a condition qu'il n'y ait pas de d\'epassement de capacit\'e, par exemple pour des entiers 32 bits, $2^{31}+2^{31}$ renverra 0. Ils sont tr\`es utilis\'es en calcul formel pour les algorithmes dits modulaires (on travaille modulo un entier assez petit). Pour travailler avec des entiers plus grands, on doit utiliser des entiers de taille plus grande, mais il faut alors programmer les op\'erations de base et d\'ecider d'un m\'ecanisme de stockage, par exemple en repr\'esentant un entier par une zone m\'emoire commencant par la taille et suivie par l'\'ecriture \`a l'aide d'entiers machines de l'entier (en base $2^{32}$). Bien entendu, plus les entiers sont grands, plus les op\'erations seront longues, par exemple l'addition de deux entiers longs de taille $N$ n\'ecessite un temps proportionnel \`a $N$, leur multiplication par l'algorithme \'el\'ementaire n\'ec\'essite un temps proportionnel \`a $N^2$ (mais il existe des algorithmes plus efficaces, par exemple Karatsuba ou FFT, cf. Knuth, The Art of Computer Programming ou la documentation de la librairie GMP). \subsection{Les rationnels.} On sait donc repr\'esenter les entiers, pour les {\bf rationnels}\index{rationnel}, il suffit de les repr\'esenter comme un couple d'entiers correspondant \`a leur \'ecriture sous forme de fraction irr\'eductible avec un d\'enominateur positif. \begin{prop} L'algorithme d'Euclide permet de calculer le PGCD (plus grand commun diviseur) de 2 entiers, \'ecrit ici en syntaxe Xcas~: \begin{verbatim} 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 } \end{verbatim} \end{prop} Preuve~: on utilise le fait qu'un nombre divise $a$ et $b$ si et seulement si il divise $r=a-bq$ et $b$. Le PGCD de $a$ et $b$ est donc le PGCD de $b$ et du reste de la division euclidienne de $a$ par $b$. Comme le reste est en valeur absolue plus petite que $|b|$, la taille des variables \verb|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\'e (ou presque, la taille des entiers longs est born\'ee parce que la taille du champ m\'emoire fixant la longueur de stockage est born\'ee) mais souvent trop lents pour les calculs num\'eriques usuels (par exemple pour calculer la valeur approch\'ee de cosinus 23 degr\'es 27 minutes). On utilise alors un autre type dont les calculs de base sont g\'er\'es par le microprocesseur (ou son coprocesseur arithm\'etique). \subsection{Les r\'eels} On se ram\`ene d'abord au cas des r\'eels positifs, en machine on garde traditionnellement un bit pour stocker le signe du r\'eel \`a repr\'esenter. \subsubsection{Virgule fixe et flottante.} La premi\`ere id\'ee qui vient naturellement serait d'utiliser un entier et de d\'eplacer la virgule d'un nombre fixe de position, ce qui revient \`a mulitplier par une puissance (n\'egative) de la base. Par exemple en base 10 avec un d\'ecalage de 4, \verb|1234.5678| serait repr\'esent\'e par \verb|12345678| et \verb|1.2345678| par \verb|12345| (on passe de l'entier au r\'eel par multiplication par $10^{-4}$. L'inconv\'enient d'une telle repr\'esentation est qu'on ne peut pas repr\'esenter des r\'eels grands ou petits, comme par exemple le nombre d'Avogadro, la constante de Planck, etc. D'o\`u l'id\'ee de ne pas fixer la position de la virgule, on parle alors de repr\'esentation \`a virgule flottante ou de nombre flottant~: on repr\'esente un nombre par deux entier, l'un appel\'e mantisse reprend les chiffres significatifs du r\'eel 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\'epare traditionnellement la mantisse de l'exposant par la lettre \verb|e|. Par exemple \verb|1234.5678| peut \^etre repr\'esent\'e par \verb|12345678e-8| (mantisse \verb|12345678|, exposant -8) mais aussi par \verb|1234567800e-10|. Naturellement, sur un ordinateur, il y a des limites pour les entiers repr\'esentant la mantisse $m$ et l'exposant $e$. Si on \'ecrit les nombres en base $b$, la mantisse $m$ s'\'ecrira avec un nombre $n$ fix\'e de chiffres (ou de bits en base 2), donc $m \in [0,b^n[$. Soit un r\'eel $x$ repr\'esent\'e par \[ x=mb^e, \quad m \in [0,b^n[ \] Si $m\in [0,b^{n-1}[$, alors on peut aussi \'ecrire $x=m' b^{e-1}$ avec $m'=mb \in [0,b^n[$, quelle \'ecriture faut-il choisir? Intuitivement, on sent qu'il vaut mieux prendre $m'$ le plus grand possible, car cela augmente le nombre de chiffres significatifs (alors que des 0 au d\'ebut de $m$ ne sont pas significatifs). Ceci est confirm\'e par le calcul de l'erreur d'arrondi pour repr\'esenter un r\'eel. En effet, si $x$ est un r\'eel non nul, il ne s'\'ecrit pas forc\'ement sous la forme $mb^e$, on doit l'arrondir, par exemple au plus proche r\'eel de la forme $mb^e$. La distance de $x$ \`a ce r\'eel est inf\'erieure ou \'egale \`a la moiti\'e de la distance entre deux flottants cons\'ecutifs, $mb^e$ et $(m+1)b^e$, donc l'erreur d'arrondi est inf\'erieure ou \'egale \`a $b^e/2$. Si on divise par $x \geq mb^e$, on obtient une erreur relative d'arrondi major\'ee par $1/(2m)$. On a donc int\'er\^et \`a prendre $m$ le plus grand possible pour minimiser cette erreur. Quitte \`a mulitplier par $b$, on peut toujours se ramener (sauf exceptions, cf. ci-dessous), \`a $m \in [b^{n-1},b^n[$, on a alors une erreur d'arrondi relative major\'ee par \[ \frac{1}{2b^{n-1}}\] On appelle {\bf flottant normalis\'e} un flottant tel que $m \in [b^{n-1},b^n[$. Pour \'ecrire un r\'eel sous forme de flottant normalis\'e, on \'ecrit le r\'eel en base $b$, et on d\'eplace la virgule pour avoir exactement $n$ chiffres non nuls avant la virgule et on arrondit (par exemple au plus proche). L'exposant est \'egal au d\'ecalage effectu\'e. 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\'e est une erreur relative inf\'erieure \`a de $2^{-53}$ ($b=2$ et $n=52+1$ pour les doubles). Exemples~: \begin{itemize} \item en base 10 avec $n=6$, pour repr\'esenter $\pi=3,14159265...$, on doit d\'ecaler la virgule de 5 positions, on obtient \verb|314159.265...| on arrondit \`a $314159$ donc on obtient \verb|314159e-5|. \item en base 2 avec $n=10$, pour repr\'esenter trois cinqui\`emes ($3/5$ en base 10, not\'e $11/101$ en base 2), on pose la division en base 2 de \verb|11| par \verb|101|, ce qui donne \begin{verbatim} 11 | 101 110 --------- -101 | 0.1001 ---- | 010 | 100 | 1000 | - 101 | ----- | 011 | \end{verbatim} on retrouve le nombre de d\'epart donc le d\'eveloppement est p\'eriodique et vaut \verb|0.1001 1001 1001 ...|. On d\'ecale le point de 10 positions, on arrondit, donc trois cinqui\`emes est repr\'esent\'e par la mantisse \verb|1001100110| et l'exposant \verb|-10|. On observe aussi sur cet exemple que $3/5$ dont l'\'ecriture en base 10 \verb|0.6| est exacte, n'a pas d'\'ecriture exacte en base 2 (de m\^eme que 1/3 n'a pas d'\'ecriture exacte en base 10). \end{itemize} Il existe une exception \`a la possibilit\'e de normaliser les flottants, lorsqu'on atteint la limite inf\'erieure de l'exposant $e$. Soit en effet $e_m$ le plus petit exposant des flottants normalis\'es et consid\'erons les flottants $x=b^{e_m}(1+1/b)$ et $y=b^{e_m}$. Ces flottants sont distincts, mais leur diff\'erence n'est plus repr\'esentable par un flottant normalis\'e. Comme on ne souhaite pas repr\'esenter $x-y$ par 0, (le test $x==y$ renvoie faux), on introduit les flottants d\'enormalis\'es, il s'agit de flottants dont l'exposant est l'exposant minimal repr\'esentable sur machine et dont la mantisse appartient \`a $[0,b^{n-1}[$. Par exemple 0 est repr\'esent\'e par un flottant d\'enormalis\'e de mantisse 0 (en fait 0 a deux repr\'sentation, une de signe positif et une de signe n\'egatif). Enfin, on utilise traditionnellement une valeur de l'exposant pour repr\'esenter les nombres plus grands que le plus grand r\'eel repr\'sentable sur machine (traditionnellement appel\'e plus ou moins infini) et les erreurs (par exemple 0./0. ou racine carr\'ee d'un nombre r\'eel n\'egatif, traditionnellement appel\'e NaN, Not a Number). {\bf Exercice}~: quels sont les nombres r\'eels repr\'esentables exactement en base 10 mais pas en base 2~? Si on écrit $1/10$ en base 2 avec 53 bits de précision, puis que l'on arrondit avec 64 bits de précision, ou si on écrit $1/10$ en base 2 avec 64 bits de précision, obtient-on la même chose~? Les ordinateurs repr\'sentent g\'en\'erallement les flottants en base 2 (cf. la section suivante pour plus de pr\'ecisions), mais cette repr\'esentation n'est pas utilis\'ee habituellement par les humains, qui pr\'ef\`erent compter en base 10. Les ordinateurs effectuent donc la conversion dans les routines d'entr\'ee-sortie. Le format standard utilis\'e pour saisir ou afficher un nombre flottant dans un logiciel scientifique est compos\'e d'un nombre \`a virgule flottante utilisant le point comme s\'eparateur d\'ecimal (et non la virgule) suivi si n\'ecessaire de la lettre \verb|e| puis de l'exposant, par exemple \verb|1.23e-5| ou \verb|0.0000123|. Dans les logiciels de calcul formel, pour distinguer un entiers repr\'esent\'es par un entier d'un entier repr\'esent\'e par un flottant on \'ecrit l'entier suivi de \verb|.0| par exemple \verb|23.0|. {\bf Remarque}~:\\ Les microprocesseurs ayant un mode BCD peuvent avoir un format de repr\'esentation des flottants en base 10, les nombres d\'ecimaux comme par exemple 0.3 peuvent \^etre repr\'esent\'es exactement. Certains logiciels, notamment maple, utilisent par d\'efaut des flottants logiciels en base 10 sur des microprocesseurs sans mode BCD, ce qui entraine une baisse de rapidit\'e importante pour les calculs num\'eriques (on peut partiellement am\'eliorer les performances en utilisant \verb|evalhf| en maple). \subsubsection{Les flottants au format {\tt double}} \index{flottant} \index{double} Cette section d\'eveloppe les notions de la section pr\'ec\'edente pour les flottants machine, utilisables dans les langage de programmation usuels, elle peut \^etre omise en première lecture. La représentation d'un double en mémoire se compose de 3 parties~: le bit\index{bit} de signe $s=\pm 1$ sur 1 bit, la mantisse\index{mantisse} $M \in [0,2^{52}[$ sur 52 bits, et l'exposant\index{exposant} $e \in [0, 2^{11}[$ sur 11 bits. Pour les nombres ``normaux'', l'exposant est en fait compris entre 1 et $2^{11}-2$, le nombre repr\'esent\'e est le rationnel \[ (1+\frac{M}{2^{52}}) 2^{e+1-2^{10}} \] Pour \'ecrire un nombre sous cette forme, il faut d'abord chercher par quel multiple de 2 il faut le diviser pour obtenir un r\'eel $r$ dans $[1,2[$, ce qui permet de d\'eterminer l'exposant $e$. Ensuite on \'ecrit la repr\'esentation en base 2 de $r-1 \in [0,1[$. Exemples~: \begin{itemize} \item -2 \\ Signe n\'egatif. Il faut diviser sa valeur absolue 2 par $2^1$ pour \^etre entre 1 et 2 dont $e+1-2^{10}=1$, l'exposant est $e=2^{10}$. On a alors $r=1$, $r-1=0$. Repr\'esentation \\ \verb|1 10000000000 00000000...0000| \item 1.5=3/2\\ Signe positif, compris entre 1 et 2 dont l'exposant v\'erifie $e+1-2^{10}=0$ soit $e=2^{10}-1=2^9+2^8+2^7+2^6+2^5+2^4+2^3+2^2+2^1+2^0$. On a $r-1=1/2=2^{-1}$. D'o\`u la repr\'esentation\\ \verb|0 01111111111 10000000...0000| \item 6.4=32/5\\ Positif. Il faut le diviser par $2^{2}$ pour avoir $8/5 \in [1,2[$ donc $e+1-2^{10}=2$ soit $e=2^{10}+1$. Ensuite $r=3/5$ qu'il faut \'ecrire en base 2 (cf. section pr\'ec\'edente), on \'ecrit donc les 52 premiers \'el\'ements du d\'eveloppement avec une r\`egle d'arrondi du dernier bit au nombre le plus proche. Ici le bit suivant le dernier \verb|1001| est un \verb|1|, on arrondit donc \`a \verb|1010|. D'o\`u la repr\'esentation\\ \verb|0 1000000001 100110011001...10011010|\\ \end{itemize} On observe que la repr\'esentation en base 2 de 6.4 a du \^etre 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\'enominateur est une puissance de 2 peuvent \^etre repr\'esent\'es exactement. Ceci entraine des r\'esultats qui peuvent surprendre comme par exemple le fait que \verb|0.3 - 3*0.1| n'est pas nul. Des représentations spéciales (avec $e=0$ ou $e=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\'esenter 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 \verb|long double| ou les types \verb|__float80| et \verb|__float128| en utilisant un drapeau de compilation du type \verb|-msse| \subsubsection{Op\'erations sur les flottants} Les opérations arithmétiques de base sur les flottants se font de la manière suivante~: \begin{itemize} \item addition et soustraction~: on détecte s'il faut additionner ou soustraire en valeur absolue en analysant les signes, on détermine l'exposant le plus grand et on décale la partie mantisse du flottant dont l'exposant est le plus petit pour se ramener à additionner deux entiers (partie mantisses correspondant au même exposant), on décale à nouveau la partie mantisse en modifiant l'exposant après l'opération pour normaliser le flottant \item multiplication~: on additionne les exposants et on multiplie les parties mantisses (vus comme des entiers), on tronque et on ajuste l'exposant si nécessaire \item division~: on soustrait les exposants et on divise les parties mantisses (division ``à virgule''), on tronque et on ajuste l'exposant si nécessaire \end{itemize} \subsubsection{Erreurs} La repr\'esentation des nombres r\'eels par des doubles pr\'esente des avantages, les op\'erations arithm\'etiques sont faites au plus vite par le microprocesseur (quoique sur les microprocesseurs plus anciens, par exemple Saturn des calculatrices HP, Z80 des calculatrices TI8x, la multiplication et la division n'est pas une opération de base du microprocesseur, elle doit être codée à partir des opérations d'addition, soustraction, décalage de bits, et/ou logique. A contrario, les coprocesseurs arithmétiques (intégrés sur les microprocesseurs de PC) proposent m\^eme 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 dans certains cas les erreurs d'arrondi). Par contre, des erreurs vont \^etre introduites, on parle de calcul approch\'e par opposition au calcul exact sur les rationnels. En effet, la repr\'esentation doit d'abord arrondir tout r\'eel qui n'est pas un rationnel dont le d\'enominateur est une puissance de 2. Ensuite chaque op\'eration va entrainer une propagation de ces erreurs et va y ajouter une erreur d'arrondi sur le r\'esultat. Enfin, l'utilisation du type double peut provoquer un d\'epassement de capacit\'e (par exemple \verb|100!*100!|). Pour diminuer ces erreurs et les risques de d\'epassement de capacit\'e, il existe des types flottants multiple pr\'ecision, qui permettent de travailler avec un nombre fix\'e \`a l'avance de d\'ecimales. Les calculs sont plus longs mais les erreurs plus faibles. Attention, il s'agit toujours de calcul approch\'e! De plus, pour des quantit\'es dont la valeur est d\'etermin\'ee de mani\`ere exp\'erimentale, la source principale de propagation d'erreurs est la pr\'ecision des quantit\'es initiales, il ne sert souvent \`a rien d'utiliser des types flottants multipr\'ecision car les erreurs dus à la représentation (double) sont négligeables devant les erreurs de mesure. \subsubsection{Erreur absolue, relative et propagation des erreurs.} On a vu pr\'ec\'edemment que pour repr\'esenter un r\'eel, on devait l'arrondir, ce qui introduit une erreur m\^eme si le r\'eel est connu exactement (par exemple 1/10). Voyons comment se propagent les {\bf erreurs}\index{erreur} 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 $|x -x_0| \leq \varepsilon_0$ et si $|y-y_0| \leq \varepsilon_1$ alors par l'inégalité triangulaire ($|a+b|\leq |a|+|b|$), on a~: \[ |(x+y)-(x_0+y_0)| \leq |x-x_0| + | y-y_0 | \leq \varepsilon_0 + \varepsilon_1 \] on dit que les erreurs {\em absolues\/} s'additionnent. \begin{defn} 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~: \[ |x-x_0| \leq \varepsilon \] \end{defn} Mais comme il faut représenter $x_0+y_0$ en machine, on doit ajouter une erreur d'arrondi, qui est proportionnelle \`a la valeur absolue de $x_0+y_0$ d'où la notion d'erreur {\em relative}~: \begin{defn} L'erreur relative est égale à l'erreur absolue divisée par la valeur absolue du nombre \[ |x-x_0| \leq \varepsilon |x_0| \] \end{defn} Remarquons au passage que les erreurs de mesure expérimentales sont pratiquement toujours des erreurs relatives. Donc lorsqu'on effectue une addition (ou une soustraction) de deux réels sur machine, on doit additionner les deux erreurs absolues sur les opérandes et ajouter une erreur d'arrondi (relative de $2^{-53}$, à titre d'exercice, on pourra vérifier que cette erreur d'arrondi est majorée par l'erreur absolue de la somme $x+y$ dès l'instant où $x$ et $y$ ont eux-même une erreur d'arrondi). Lorsqu'on effectue une multiplication de deux nombres $x,y$ dont les représentants $x_0,y_0$ sont non nuls, on a \[ \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^{-53}$ sur $x_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 $x$ et $y$ non nuls est de $2^{-53}$, alors l'erreur relative sur $xy$ sera de \[ 2^{-53} + 2^{-53} + 2^{-106} + 2^{-53} \approx 3 \times 2^{-53} \] Lorsque l'erreur relative sur les données est grande devant $2^{-53}$, l'erreur relative d'arrondi final est négligeable, on peut alors dire que les erreurs relatives s'additionnent pour un produit (c'est aussi vrai pour un quotient: exercice!). Par contre, si on additionne deux nombres dont le représentant de la somme est proche de 0, la somme des erreurs absolues peut devenir non négligeable par rapport à la somme des représentants, entrainant une erreur relative très grande. Par exemple si $x$ est représenté par $x_0=1+2^{-52}$ avec une erreur d'arrondi de $2^{-53}$ et $y$ par $y_0=-1$ avec la même erreur d'arrondi, l'addition de $x$ et $y$ renvoie $2^{-52}$ avec une erreur absolue de $2 * 2^{-53}$ (ici il n'y a pas d'arrondi lorsqu'on fait la somme). C'est une erreur relative de $1$ (qui domine largement l'erreur d'arrondi) ce qui signifie que dans la mantisse, seul le premier bit sur les 52 a un sens, la perte de précision est très grande. Une autre conséquence importante est que l'addition de réels sur machine n'est pas une opération associative, par exemple \[ (2.0^{-53}+2.0^{-53})+1.0 \rightarrow 1+2^{-52} \] alors que \[ (2.0^{-53}+1.0)+2.0^{-53} \rightarrow 1 \] 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\'erique d'une d\'eriv\'ee de fonction, il vaut mieux calculer $(f(x+h)-f(x-h))/(2h)$ que $(f(x+h)-f(x))/h$. Attention \`a ne pas prendre $h$ trop petit, sinon $x+h=x$. 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. Il est d'ailleurs souvent trop difficile de calculer la majoration rigoureuse de l'erreur pour des calculs complexes. 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. On peut aussi faire varier l\'eg\`erement les donn\'ees et observer la sensibilit\'e du r\'esultat. 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 (par exemple la bibliothèque C MPFI). \subsection{Types composés.} Après les nombres réels, on passe aux {\bf nombres complexes}\index{complexe}~: on utilise un couple (partie réelle, imaginaire) de fractions (exacts) ou de flottants et les r\`egles habituelles sur les complexes. Après les nombres, l'objet le plus utilisé dans les systèmes de calcul formel est probablement le {\bf polynôme}\index{polynome}, 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~: \begin{itemize} \item les polynômes à 1 variable, représentation dense, on stocke la liste des coefficients du polynôme par ordre croissant ou décroissant \item les polynômes à 1 variable, représentation creuse, on stocke des paires coefficients, degré pour les coefficients non nuls \item les polynômes à plusieurs variables, représenté récursivement de manière dense ou creuse (i.e. $P(x_1,...,x_n)$ vu comme polynôme en $x_n$ à coefficients polynômes dépendant des variables $x_1,...,x_{n-1}$), ce sont des cas particuliers des 2 cas précédents \item les polynômes à plusieurs variables distribués, on stocke des monômes, qui sont des paires coefficient, liste d'entiers, la liste représentant les exposant des variables dans le monôme. \item la représentation symbolique (par exemple $x y^2-5x+y^3$) beaucoup plus difficile à manipuler directement \end{itemize} 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\^eme algorithme que pour les entiers mais avec la division euclidienne des polyn\^omes, il existe des algorithmes plus efficaces, cf. le chapitre sur les polyn\^omes) Lien avec la repr\'esentation en base $z$ (TD). Les polyn\^omes peuvent servir \`a repr\'esenter des nombres non rationnels de mani\`ere exacte, par exemple les nombres alg\'ebriques (qui sont solutions d'une \'equation polynomiale). Les {\bf symboles}\index{symbole} ou noms de variable d\'esignent par exemple le nom d'une inconnue dans un polyn\^ome, ils sont repr\'esent\'es par une chaine de caract\'ere et peuvent \^etre affect\'es \`a 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 {\bf expressions}\index{expression}, par exemple \verb|sin(x)+2*x^2|, elles peuvent \^etre repr\'esent\'ees 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 (\verb|''|) ou implicitement (par exemple l'affectation n'évalue pas le symbole qu'on va affecter). Les {\bf fonctions}\index{fonction} ne doivent pas \^etre confondues avec les expressions, elles associent \`a leurs arguments une expression. Par exemple \verb|sin| est une fonction, alors que \verb|sin(x)| est une expression. Les conteneurs contiennent plusieurs objets et permettent d'associer \`a un indice un objet. Il en existe de plusieurs types, par exemple les {\bf listes}\index{liste} et les {\bf s\'equences}\index{sequence} dont l'indice est un entier compris entre 1 (ou 0) et la taille (-1), les {\bf tables} dont l'indice est plus g\'en\'eral, et les tableaux (utilisés pour les {\bf vecteurs, matrices}\index{vecteur}\index{matrice}) 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 \verb|[]|. Il faut prendre garde au fait qu'en général affecter par exemple \verb|l[1]:=3;| à une variable libre \verb|l| crée une table et non une liste. Remarque: certains logiciels acc\'edent \`a certains types de conteneurs uniquement par r\'ef\'erence (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 \'el\'ement du conteneur modifie alors toutes les copies pointant sur ce conteneur. Cette m\'ethode est plus efficace mais peut \^etre surprenante. %Idées TD: Exercices sur les nombres en base 2, %sur les changements de base, sur les erreurs. %Ecriture en base $z$ avec reste symétrique. %L'addition n'est pas associative dans les flottants : comment %sommer des séries de terme général décroissant pour avoir la meilleure %précision possible. %Changement d'origine pour un polynome avec appel de Taylor. \section{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. \subsection{Le point fixe} Soit $f$ une fonction continue sur un intervalle $I=[a,b]$ de $\R$, et à valeurs dans $I$ (attention à bien choisir $I$ pour que l'image de $I$ par $f$ reste dans $I$). On s'intéresse à la suite \begin{equation} \label{eq:u_n_fixe} u_{n+1}=f(u_n), \quad u_0 \in I \end{equation} Supposons que $u_n$ converge vers une limite $l \in I$ lorsque $n \rightarrow +\infty$, alors la limite doit vérifier \[ f(l)=l \] puisque $f$ est continue. On dit que $l$ est un point fixe de $f$. Ceci amène à l'idée d'utiliser ces suites pour résoudre numériquement l'équation $f(x)=x$. Nous allons donner un théorème permettant d'assurer que la suite (\ref{eq:u_n_fixe}) converge, et que la limite est l'unique solution de $f(l)=l$ sur $I$. \begin{defn} On dit que $f$ est contractante\index{contractante} de rapport $k<1$ sur $I$ si \[ \forall x,y \in I, \quad |f(y)-f(x)| \leq k |y-x| \] \end{defn} En pratique, les fonctions $f$ que l'on considèrera seront continument dérivables, donc d'après le théorème des accroissements finis \[ f(y)-f(x)=f'(\theta) (y-x), \quad \theta \in [x,y] \] ainsi pour vérifier que $f$ est contractante, on étudie la valeur absolue de $f'$ sur $I$, il suffit de montrer que cette valeur absolue est strictement inférieure à un réel $k<1$ pour conclure (il faut donc chercher le maximum de $|f'|$ sur $I$. Attention, il s'agit du maximum de $|f'|$ et pas du maximum de $f'$, ce qui revient à chercher le maximum de $f'$ et de $-f'$). On a alors le \begin{theorem} (du point fixe)\index{point fixe}\\ si $f$ est contractante de $I=[a,b]$ dans $I$ de rapport $k$ alors la suite (\ref{eq:u_n_fixe}) converge vers l'unique solution de $f(l)=l$ dans $I$. On a de plus les encadrements~: \begin{equation} \label{eq:u_n_l} |u_n-l| \leq k^n |b-a|, \quad |u_n -l | \leq \frac{|u_{n+1}-u_n|}{1-k} \end{equation} \end{theorem} Démonstration~: Tout d'abord si $f$ est contractante, on montre à partir de la définition de la continuité que $f$ est continue. Soit $g(x)=f(x)-x$, alors $g$ est continue, positive en $a$ et négative en $b$, il existe donc $l\in[a,b]$ tel que $g(l)=0$ (théorème des valeurs intermédiaires). Soit $u_n$ une suite définie par (\ref{eq:u_n_fixe}). On a alors pour tout $n$ \[ |u_{n+1}-l|=|f(u_n)-f(l)| \leq k |u_n-l| \] Donc par une récurrence évidente~: \[ |u_n-l| \leq k^n |u_0-l| \] ce qui entraine d'ailleurs que $|u_n-l| \leq k^n |a-b|$. Comme $k \in [0,1[ $, la suite géométrique $k^n$ converge vers 0 lorsque $n$ tend vers l'infini, donc $u_n$ tend vers $l$. Notons que $l$ est unique car si $l'$ est une autre solution alors $|l-l'|=|f(l)-f(l')| \leq k|l-l'|$ donc $(1-k)|l-l'| \leq 0$, or $1-k>0$ et $|l-l'| \geq 0$ donc $|l-l'|$ doit \^etre nul. Passons à la preuve de la majoration (\ref{eq:u_n_l}) qui est importante en pratique car elle donne un test d'arrêt de calcul des termes de la suite récurrente, on écrit pour $m>0$~: \[ 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_n-l| \leq \sum_{j=0}^{m-1} |u_{n+j}-u_{n+j+1}| + |u_m-l| \] puis on applique le fait que $f$ est contractante de rapport $k$ \[ |u_n-l| \leq \sum_{j=0}^{m-1} k^j |u_{n}-u_{n+1}| + |u_m-l| \] soit \[ |u_n-l| \leq \frac{1-k^m}{1-k} |u_{n}-u_{n+1}| + |u_m-l| \] On fait alors tendre $m$ vers l'infini d'où le résultat. Exemples~: Cherchons une valeur approchée de $\sqrt{2}$ par cette méthode. Il faut d'abord trouver une fonction $f$ dont $\sqrt{2}$ est un point fixe, par exemple \[ f(x)=\frac{x+2}{x+1}\] On vérifie que $f(\sqrt{2})=\sqrt{2})$, puis que $f'=-1/(x+1)^2$ donc $f$ décroit. On va voir si les hypothèses du théorème du point fixe s'appliquent sur par exemple $[1,2]$. Comme $f$ est décroissante $f([1,2])=[f(2),f(1)]=[4/3,3/2]$ qui est bien inclus dans $[1,2] $. De plus $f'$ est comprise entre $-1/(1+1)^2=-1/4$ et $-1/(2+1)^2=-1/9$ donc $|f'|<1/4$, $f$ est contractante de rapport $1/4$. On peut donc itérer la suite à partir par exemple de $u_0=1$ et on va converger vers $\sqrt{2}$ (en s'en rapprochant à chaque cran d'un rapport inférieur à $1/4$). Considérons l'équation en $x$ \[ x- e \sin(x) =t, \quad e \in [0,1[ \] c'est l'équation du temps utilisée en astronomie pour trouver la position d'une planète sur son orbite elliptique ($e$ étant l'excentricité de l'ellipse). Il n'y a pas de formule exacte permettant de calculer $x$ en fonction de $t$. Si on a une valeur numérique pour $t$, on peut trouver une valeur numérique approchée de $x$ par la méthode du point fixe, en réécrivant l'équation sous la forme \[ f(x)=t+e\sin(x) = x \] On observe que $f$ envoie $\R$ dans $[t-e,t+e]$ donc on peut prendre $I=[t-e,t+e]$, de plus $|f'|\leq e <1$, $f$ est contractante de rapport $e \in [0,1[$, le théorème s'applique, il suffit de prendre une valeur initiale dans $[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 $x$ à $10^{-6}$ près, il suffira que la différence entre deux termes successifs de la suite $u_n$ vérifie \[ |u_{n+1}-u_n| \leq 10^{-6} (1-e) \] on aura alors bien~: \[ |u_n-x| \leq \frac{|u_{n+1}-u_n|}{1-e} \leq 10^{-6} \] Cette méthode n'est pas toujours optimale, car la vitesse de convergence vers la limite $l$ est dite ``linéaire'', c'est-\`a-dire que le temps de calcul pour avoir $n$ décimales est proportionnel à $n$ (ou encore il faut effectuer un nombre d'itérations proportionnel à $n$, chaque itération faisant gagner en précision de l'ordre du rapport $k$ de contractance). En effet, supposons que $f'$ est continue en $l$ et que $00$ tel que la suite définie par \[ u_{n+1} = u_n -\frac{f(u_n)}{f'(u_n)}, \quad |u_0-r| \leq \varepsilon \] converge vers $r$. Si on a $|f'{'}| \leq M$ et $|1/f'| \leq m$ sur un intervalle $[r-\eta,r+\eta]$ contenu dans $I$, alors on peut prendre tout r\'eel $\varepsilon>0$ tel que $\varepsilon < 2/(mM)$ et $\varepsilon \leq \eta$. \end{theorem} {\bf Démonstration}~: on a \[ 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 $f$ en $u_n$ à l'ordre 2, on obtient pour un réel $\theta$ situé entre $r$ et $u_n$~: \[ 0 = f(r)=f(u_n)+(r-u_n) f'(u_n) + (r-u_n)^2 \frac{f'{'}(\theta)}{2} \] donc~: \[ (u_n-r)f'(u_n)-f(u_n)= (u_n-r)^2 \frac{f'{'}(\theta)}{2} \] d'o\`u~: \[ |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-\varepsilon,r+\varepsilon]$ contenant strictement $r$ et tel que $|f'{'}|x_0$, donc $g$ est décroissante pour $xx_0$. On conclut alors que $g \geq 0$ puisque $g(x_0)=0$. Donc $f$ est bien au-dessus de sa tangente. On arrive au deuxième théorème sur la méthode de Newton \begin{theorem} \index{Newton} Si $f(r)=0, f'(r)>0$ et si $f'{'} \geq 0$ sur $[r,b]$ alors pour tout $u_0 \in [r,b]$ la suite de la méthode de Newton \[ u_{n+1} = u_n -\frac{f(u_n)}{f'(u_n)}, \] est définie, décroissante, minorée par $r$ et converge vers $r$. De plus \[ 0 \leq u_n -r \leq \frac{f(u_n)}{f'(r)} \] \end{theorem} {\bf Démonstration}~:\\ On a $f'{'} \geq 0$ donc si $f'(r)>0$ alors $f'>0$ sur $[r,b]$, $f$ est donc strictement croissante sur $[r,b]$ on en déduit que $f>0$ sur $]r,b]$ donc $u_{n+1} \leq u_n$. Comme la courbe représentative de $f$ est au-dessus de la tangente, on a $u_{n+1} \geq r$ (car $u_{n+1}$ est l'abscisse du point d'intersection de la tangente avec l'axe des $x$). La suite $u_n$ est donc décroissante minorée par $r$, donc convergente vers une limite $l\geq r$. \`A la limite, on a \[ l=l-\frac{f(l)}{f'(l)} \Rightarrow f(l)=0 \] donc $l=r$ car $f>0$ sur $]r,b]$. Comme $(u_n)$ est d\'ecroissante, on a bien $0 \leq u_n -r$, pour montrer l'autre in\'egalit\'e, on applique le th\'eor\`eme des accroissements finis, il existe $\theta \in [r,u_n]$ tel que \[ f(u_n)-f(r)=(u_n-r)f'(\theta) \] comme $f(r)=0$, on a \[ u_n-r = \frac{f(u_n)}{f'(\theta)} \] et la deuxi\`eme in\'egalit\'e du th\'eor\`eme en d\'ecoule parce que $f'$ est croissante. {\bf Variantes}~:\\ Il existe des variantes, par exemple si $f'(r)<0$ et $f'{'} \geq 0$ sur $[a,r]$. Si $f'{'} \leq 0$, on considère $g=-f$. {\bf Application}~:\\ On peut calculer la valeur approch\'ee de la racine $k$-i\`eme d'un r\'eel $a>0$ en appliquant ce deuxi\`eme th\'eor\`eme. En effet si $a>0$, alors $x^k-a$ est 2 fois continument d\'erivable et de d\'eriv\'ee premi\`ere $kx^{k-1} $ et seconde $k(k-1)x^{k-2}$ strictement positives sur $\R^{+\ast}$ (car $k \geq 2$). Il suffit donc de prendre une valeur de d\'epart $u_0$ plus grande que la racine $k$-i\`eme, par exemple $1+a/k$ (en effet $(1+a/k)^k \geq 1+k a/k=1+a$). En appliquant l'in\'egalit\'e du th\'eor\`eme, on a~: \[ 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\'ee de $\sqrt[k]{a}$ \`a $\varepsilon$ pr\`es, on peut donc choisir comme test d'arr\^et \[ u_n^k -a \leq \frac{ka}{1+\frac{a}{k}} \varepsilon \] Par exemple pour $\sqrt{2}$, le test d'arr\^et serait $u_n^2-2 \leq 2 \varepsilon$. %% Idées TD: différentes suites itératives approchant $\sqrt{3}$, %% par exemple trouver (a*x+b)/(c*x+d) ayant sqrt(3) comme point fixe %% vitesse de convergence théorique. %% Recherche de suite itérative pour les racines $n$-ième. %% Idées TP: observation pratique de la vitesse de convergence, %% en précision double ou illimitée. Exemple de cas où Newton marche %% et ne marche pas, utilisation du préfacteur. \section{Développement de Taylor, s\'eries enti\`eres, fonctions usuelles} Résumé: Séries entières. Calcul des fonctions transcendantes usuelles. Soit $f$ une fonction indéfiniment dérivable sur un intervalle $I$ de $\R$ et $x_0 \in I$. On peut alors effectuer le développement de Taylor de $f$ en $x_0$ à l'ordre $n$\index{Taylor} \[ 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)$ converge lorsque $n$ tend vers l'infini, si la limite est égale à $f(x)$ et si on peut facilement majorer la différence entre $f(x)$ et $T_n(f)(x)$. Si c'est le cas, on pourra utiliser $T_n(f)(x)$ comme valeur approchée de $f(x)$. On peut parfois répondre à ces questions simultanément en regardant le développement de Taylor de $f$ avec reste~: il existe $\theta$ compris entre $x_0$ et $x$ tel que \[ R_n(x) := f(x)- T_n(f)(x) = (x-x_0)^{n+1}\frac{f^{[n+1]}(\theta)}{(n+1)!} \] C'est le cas pour la fonction exponentielle que nous allons détailler, ainsi que les fonctions sinus et cosinus. \subsection{La fonction exponentielle}\index{exp} Soit $f(x)=\exp(x)$ et $x_0=0$, la dérivée $n$-ième de $f$ est $\exp(x)$, donc $R_n(x)=\exp(\theta)x^{n+1}/(n+1)!$ avec $\theta$ compris entre 0 et $x$, ainsi si $x$ est positif $|R_n(x)| \leq e^x x^{n+1}/(n+1)!$ et si $x$ est négatif, $|R_n(x)| \leq x^{n+1}/(n+1)!$. Dans les deux cas, la limite de $R_n$ est 0 lorsque $n$ tend vers l'infini, car pour $n \geq 2x$, on a \[ \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 $x$ réel \[ 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!} \] Comment en déduire une valeur approchée de $e^x$? Il suffira d'arr\^eter la sommation lorsque $R:=x^{n+1}/(n+1)!$ si $x<0$ ou lorsque $R:=e^x x^{n+1}/(n+1)!$ si $x>0$ est inférieur à l'erreur absolue souhaitée, le plus tôt étant le mieux pour des raisons d'efficacité et pour éviter l'accumulation d'erreurs d'arrondi. Si on veut connaitre $e^x$ à une erreur relative $\varepsilon$ donnée (par exemple $\varepsilon=2^{-53}$ pour stocker le résultat dans un double) il suffit que $R/e^x < \varepsilon$, donc si $x$ est positif, il suffit que $x^{n+1}/(n+1)!<\varepsilon$, on peut donc arrêter la sommation lorsque le terme suivant est plus petit que $\varepsilon$. On observe que plus $x$ est grand, plus $n$ devra être grand pour réaliser le test d'arrêt, ce qui est facheux pour le temps de calcul. De plus, le résultat final peut être petit alors que les termes intermédiaires calculés dans la somme peuvent être grands, ce qui provoque une perte de précision relative, par exemple si on veut calculer $e^{-10}$ ou plus géné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^{-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^x$ (ne change pas l'erreur relative). Pour les grands réels, on peut utiliser $e^{2x}=(e^x)^2$ (multiplie par 2 l'erreur relative). On peut aussi, si on connait une valeur approchée de $\ln(2)$, effectuer la division euclidienne de $x$ par $\ln(2)$ avec reste symétrique~: \[ x = a \ln(2) + r, \quad a \in \Z, |r| \leq \frac{ln(2)}{2} \] puis si $r$ est positif, on somme la série de $T(f)(r)$, si $r$ est négatif, on calcule $T(f)(-r)$ et on inverse, on applique alors~: \[ e^x = 2^a e^r \] Il faut toutefois noter que $\ln(2)$ n'étant pas connu exactement, on commet une erreur d'arrondi absolu sur $r$ d'ordre $a \eta$, où $\eta$ est l'erreur relative sur $\ln(2)$, il faut donc ajouter une erreur d'arrondi relative de $x/\ln(2) \eta$ qui peut devenir grande si $x$ est grand. Puis il faut ajouter la somme des erreurs d'arrondi due au calcul de $e^r$, que l'on peut minimiser en utilisant la méthode de Horner pour évaluer $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)$ avec cette précision, le calcul de $e^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^{-53}$). Notons que en général $x$ lui-même a déjà été arrondi ou n'est connu qu'avec une précision relative. Or si $x>0$ est connu avec une erreur relative de $\varepsilon$ (donc une erreur absolue de $\varepsilon |x|$, alors \[ e^{x+\varepsilon |x|}= e^x e^{\varepsilon |x|} \] donc on ne peut pas espérer mieux qu'une erreur relative de $e^{\varepsilon |x|}-1$ sur l'exponentielle de $x$. Si $\varepsilon x$ est petit cette erreur relative (impossible à éviter, quel que soit l'algorithme utilisé pour calculer l'exponentielle) est d'ordre $\varepsilon |x|$. Si $\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 $x$ est supérieur à 709, l'exponentielle renvoie infini). Exercice~: refaire les mêmes calculs pour les fonction sinus ou cosinus. On utilise par exemple $\sin(x+\pi)=-sin(x)$, $\sin(-x)=-\sin(x)$, $\sin(x)=\cos(\pi/2-x)$ pour se ramener au calcul de $\sin(x)$ ou de $\cos(x)$ sur $[0,\pi/4]$.\index{sin}\index{cos} \[ \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 $n$-ième d'une fonction (par exemple avec $\tan(x)$), et encore plus de la majorer. D'où l'intérêt de développer une théorie des fonctions qui sont égales à leur développement de Taylor à l'infini d'une part, et d'avoir d'autres méthodes pour majorer le reste, nous présentons ici le cas des séries alternées. \subsection{Séries entières.}\index{serie entiere} Les séries de type prendre la limite lorsque $n$ tend vers l'infini du développement de Taylor en x=0 sont de la forme \[ \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 à $\sum_{n=0}^\infty a_n x^n$ lorsque $a_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_0$ tel que $|a_n x_0^n|$ est born\'e (ce sera le cas en particulier si la s\'erie converge en $x_0$), alors \[ |a_n x^n| = |a_n x_0^n| |\frac{x}{x_0}|^n \leq M |\frac{x}{x_0}|^n \] la s\'erie converge donc en $x$ si $|x|<|x_0|$ et on peut majorer le reste de la s\'erie au rang $n$ par \[ |R_n| \leq M \frac{ |\frac{x}{x_0}|^n} {1-|\frac{x}{x_0}|} \] la vitesse de convergence est donc du m\^eme type que pour le th\'eor\`eme du point fixe (le nombre de termes \`a calculer pour trouver une valeur approch\'ee avec $k$ d\'ecimales d\'epend lin\'eairement $k$, les constantes sont d'autant plus grandes que $|x|$ est grand). \begin{theorem} S'il existe un rang $n_0$, un r\'eel $M>0$ et un complexe $x_0$ tels que pour $n>n_0$, on ait~: \[ |a_n x_0|^n \leq M\] alors la s\'erie converge pour $|x|<|x_0|$ et pour $n\geq n_0$, on a~: \begin{equation} \label{eq:maj_serie_entiere} |R_n| \leq M \frac{ |\frac{x}{x_0}|^n} {1-|\frac{x}{x_0}|} \end{equation} \end{theorem} On en d\'eduit qu'il existe un réel positif $R\geq 0$ éventuellement égal à $+\infty$ tel que la s\'erie converge (la limite de la somme jusqu'\`a l'infini existe) lorsque $|x|R$, ce réel est appelé {\bf 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 $\sum x^n$ (car elle diverge si $|x|>1$ et converge si $|x|<1$). On ne peut pas dire ce qui se passe g\'en\'eriquement lorsqu'on est \`a la limite, c'est-\`a-dire lorsque $|x|=R$ (si $R\neq +\infty$). Mais cela n'a en fait pas trop d'importance en pratique car m\^eme si la s\'erie converge, elle converge souvent trop lentement pour donner de bonnes approximations. En fait, la vitesse de convergence d'une s\'erie enti\`ere de rayon $R\neq +\infty$ est en gros la m\^eme que celle d'une s\'erie g\'eom\'etrique de raison $|x|/R$. Lorsque 2 séries ont un rayon de convergence non nul, alors on peut effectuer leur somme, leur produit comme des polynômes et la série somme/produit a un rayon de convergence au moins égal au plus petit des 2 rayons de convergence des arguments. On peut inverser une série entière non nulle en 0 en appliquant \[ (1+x)^{-1} = 1-x+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 $g$ et $f$ en $g\circ f$ (avec les r\`egles de calcul de composition des polyn\^omes) si $f(0)=0$. On peut enfin dériver et intégrer une série entière terme à terme dans son rayon de convergence. On dit qu'une fonction est développable en série entière en 0 si elle est égale à son développement de Taylor en 0 sommé jusqu'en l'infini dans un disque de centre 0 et de rayon non nul. Les fonctions exponentielle, sinus, cosinus sont donc développables en série entière en 0. La fonction tangente également car le dénominateur cosinus est non nul en 0, mais son rayon de convergence n'est pas l'infini et le calcul des $a_n$ est assez complexe. La fonction $(1+x)^\alpha$ est développable en séries entières pour tout $\alpha \in \R$ avec un rayon de convergence 1 (ou l'infini pour $\alpha $ entier positif). \[ (1+x)^\alpha = 1 + \alpha x + \frac{\alpha (\alpha-1)}{2!} x^2 + ... + \frac{\alpha (\alpha-1) ... (\alpha -n +1)}{n!} x^n + ...\] Pour $\alpha=-1$, c'est la série géométrique de raison $-x$, en effet si $|x|<1$~: \[ \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 à $x$, on obtient que $\ln(1+x)$ est développable en série entière en 0 de rayon de convergence 1 et \[ \ln(1+x) = \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)$ en iintégrant celui de $1/(1+x^2)$, de même pour $\arccos(x)$ et $\arcsin(x)$ en intégrant celui de $(1-x^2)^{-1/2}$.\index{atan} \[ \arctan(x)=\sum_{n=0}^\infty (-1)^{n} \frac{x^{2n+1}}{2n+1}, \] On peut donc calculer $\ln$, $\arctan$, ... par ces formules, mais il faut répondre à la question où arrête-t-on la somme pour obtenir une précision donnée? Dans le cas de $\ln(1+x)$, on pourrait répondre comme avec l'exponentielle en majorant la dérivée $n+1$-ième, mais ce n'est plus faisable pour $\arctan, \arcsin, \arccos$. On va donner un autre critère qui ne nécessite pas de calculer cette dérivée mais utilise l'alternance des signes dans la somme. \subsection{Série alternée}\index{serie alternee} \begin{theorem} Soit $S_n= \sum_{k=0}^n (-1)^k u_k$ la somme jusqu'au rang $n$ d'une s\'erie de r\'eels tels que la suite des $u_k$ d\'ecroit \`a partir d'un rang $n_0$ et tend vers 0 lorsque $k\rightarrow +\infty$. Alors $S_n$ converge vers une limite $S$. Si $n\geq n_0$, la limite est comprise entre deux sommes partielles succesives $S_n$ et $S_{n+1}$ et le reste est major\'e par la valeur absolue du premier terme non somm\'e~: \[ |R_n| \leq |u_{n+1}|\] \end{theorem} D\'emonstration~:\\ on montre que les suites $v_n=S_{2n}$ et $w_n=S_{2n+1}$ sont adjacentes. On a \[ 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_n$ est d\'ecroissante, de m\^eme $w_n$ est croissante, et $v_n-w_n=u_{2n+1}$ est positif et tend vers 0. On en d\'eduit que $v_n$ et $w_n$ convergent vers la m\^eme limite $S$ telle que $v_n>S>w_n$ et les in\'egalit\'es du th\'eor\`eme s'en d\'eduisent. {\bf Remarque}\\ lorsqu'on utilise une suite altern\'ee pour trouver une valeur approch\'ee, il faut que $u_n$ tende assez vite vers 0, sinon il y aura perte de pr\'ecision sur la mantisse lorsqu'on effectuera $u_{2n}-u_{2n+1}$. On sommera aussi les termes par ordre d\'ecroissant pour diminuer les erreurs d'arrondi. \subsection{La fonction logarithme}\index{ln} Si nous voulons calculer $\ln(1+x)$ pour $x \in [0,1[$ avec une pr\'ecision $\varepsilon$, il suffit de calculer \[ \sum_{k=0}^n (-1)^k \frac{x^{k+1}}{k+1}\] pour $n$ tel que la valeur absolue du terme suivant soit plus petit que $\varepsilon$~: \[ n \mbox{ tel que } \frac{x^{n+1}}{n+1} < \varepsilon \] en effet, les signes sont altern\'es et la suite $\frac{x^{k+1}}{k+1}$ d\'ecroit vers 0. Si la suite d\'ecroit lentement vers 0, cette m\'ethode est mauvaise num\'eriquement et en temps de calcul car il y a presque compensation entre termes successifs donc perte de pr\'ecision sur la mantisse et il y a beaucoup de termes \`a calculer. C'est le cas pour le logarithme, si $x$ est voisin de 1, il faut calculer $n$ termes pour avoir une pr\'ecision en $1/n$, par exemple 1 million de termes pour avoir une pr\'ecision de $1e-6$ (sans tenir compte des erreurs d'arrondi). Si $x$ est proche de $1/2$ il faut de l'ordre de $-\ln(\varepsilon)/\ln(2)$ termes ce qui est mieux, mais encore relativement grand (par exemple 50 termes environ pour une pr\'ecision en $1e-16$, 13 termes pour $1e-4$). On a donc int\'er\^et \`a se ramener si possible \`a calculer la fonction en un $x$ o\`u la convergence est plus rapide (donc $|x|$ le plus petit possible). Par exemple pour le calcul de $\ln(1+x)$ on peut~: \begin{itemize} \item utiliser la racine carr\'ee \[ ln(1+x)= 2 ln(\sqrt{1+x})\] on observe que~: \[ X=\sqrt{1+x}-1 = \frac{x}{1+\sqrt{1+x}} \leq \frac{x}{2} \] il faut toutefois faire attention \`a la perte de pr\'ecision sur $X$ par rapport \`a $x$ lorsque $x$ est petit. \item utiliser l'inverse \[ \ln(1+x)=-\ln(1/(1+x))=-\ln(1 + \frac{-x}{1+x}) \] lorsque $x$ est proche de 1, $-x/(1+x)$ est proche de $-x/2$, on a presque divis\'e par 2. Attention toutefois, on se retrouve alors avec une s\'erie non altern\'ee, mais on peut utiliser (\ref{eq:maj_serie_entiere}) pour majorer le reste dans ce cas. \item trouver une valeur approch\'ee $y_0$ de $\ln(1+x)$ \`a une pr\'ecision faible, par exemple $1e-4$, et utiliser la m\'ethode de Newton pour am\'eliorer la pr\'ecision. Soit en effet $y=\ln(1+x)$, alors $e^y=1+x$, on pose $f(y)=e^y-(1+x)$, on utilise la suite it\'erative \[ y_{n+1} = y_n - \frac{e^{y_n}-(1+x)}{e^{y_n}} \] Comme $y_0$ est proche \`a $1e-4$ de $y$, on peut esp\'erer avoir une valeur approch\'ee de $y$ \`a $1e-16$ en 2 it\'erations. Notez que $y$ est proche de $0$, on est dans un domaine o\`u le calcul de $e^y$ est rapide et pr\'ecis et de plus la m\'ethode de Newton ``corrige'' les erreurs interm\'ediaires. \end{itemize} Nous sommes donc en mesure de calculer pr\'ecis\'ement le logarithme $\ln(1+x)$ pour disons $|x|<1/2$. Pour calculer $ln$ sur $\R^+$, on se ram\`ene \`a $[1,2]$ en utilisant l'\'ecriture mantisse-exposant, puis si $x\in[3/2,2]$ on peut en prendre la racine carr\'ee pour se retrouver dans l'intervalle souhait\'e. On peut aussi effectuer une division par $\sqrt{2}$. Remarquons que si $x$ est connu \`a une erreur relative $\varepsilon$ pr\`es, comme \[ \ln(x(1 \pm \varepsilon))=\ln(x) + \ln(1 \pm \varepsilon) \] $\ln(x)$ est connu \`a une erreur absolue de $|\ln(1 \pm \varepsilon)| \approx \varepsilon$. Si $\ln(x)$ est proche de 0, on a une grande perte de pr\'ecision relative. Finalement, nous savons calculer $\ln$ et $\exp$ sous r\'eserve d'avoir dans une table la valeur de $\ln(2)$. Pour calculer $\ln(2)$ pr\'ecis\'ement, on peut utiliser \[ \ln(2)=-ln(1/2)=-ln(1-1/2) \] et le d\'eveloppement en s\'erie calcul\'e en mode exact avec des fractions \`a un ordre suffisant, on majore le reste en utilisant que le terme g\'en\'eral de la s\'erie $\ln(1+x)$ est born\'e par $M=1$ en $x=1$, donc d'apr\`es (\ref{eq:maj_serie_entiere})~: \[ |R_n| \leq \frac{1}{2^n}\] (on peut m\^eme obtenir $1/(n2^n)$ car on a besoin de $M$ uniquement pour les termes d'ordre plus grand que $n$, on peut donc prendre $M=1/n$). Par exemple, pour avoir $\ln(2)$ avec une mantisse de 80 bits, on effectue une fois pour toutes avec un logiciel de calcul formel~:\\ \verb|a:=sum((1/2)^k/k,k=1..80)||\\ puis la division en base 2 avec 81 bits de pr\'ecision \verb|iquo(numer(a)*2^81,denom(a))| Exercice~: pour les fonctions trigonom\'etriques, il faut une m\'ethode de calcul de $\pi$. On peut par exemple faire le calcul de $16 \arctan(1/5)-4\arctan(1/239)$ en utilisant le d\'eveloppement de la fonction $\arctan$ \`a un ordre suffisant. \subsection{Autres applications} On peut calculer certaines int\'egrales de la m\^eme mani\`ere, par exemple \[ \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). \subsubsection{Exemple~: la fonction d'erreur (error fonction, {\tt erf})} Cette fonction est définie à une constante multiplicative près par~: \[ 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)= \sum_{n=0}^{+\infty} (-1)^n \frac{x^{2n}}{n! (2n+1)}\] Ce développement converge très rapidement pour $|x|\leq 1$. Par contre, pour $|x|$ grand, il faut calculer beaucoup de termes avant que le reste soit suffisamment petit pour être négligeable, et certains termes intermédiaires sont grands, ce qui provoque une perte de précision qui peut rendre le résultat calculé complètement faux. Contrairement à la fonction exponentielle, il n'y a pas de possibilité de réduire l'argument à une plage où la série converge vite. Il faut donc \begin{itemize} \item soit utiliser des flottants multiprécision, avec une précision augmentée de la quantité nécessaire pour avoir un résultat fiable \item soit, pour les grandes valeurs de $x$, utiliser un développement asymptotique (en puissances de $1/x$) de \[ \int_x^{+\infty} e^{-t^2} \ dt \] ainsi que \[ \int_0^{+\infty} e^{-t^2} \ dt =\frac{\sqrt{\pi }}{2} \] Le développement asymptotique s'obtient par exemple en changeant de variable $u=t^2$ et en effectuant des intégrations par parties répétées en intégrant $e^{-u}$ et en dérivant $u^{-1/2}$ et ses dérivées successives. Ce type de développement asymptotique a la propriété inverse du développement en 0: les termes successifs commencent par décroitre avant de croitre et de tendre vers l'infini. Il faut donc arrêter le développement à un rang donné (dépendant de $x$) et il est impossible d'obtenir une précision meilleure pour cette valeur de $x$ par un développement asymptotique (on parle parfois de développement des astronomes). \end{itemize} {\bf Exercice~}: donner une valeur approchée de $f(1)$ à $1e-16$ près. Combien de termes faut-il calculer dans la somme pour trouver une valeur approchée de $f(7)$ à $1e-16$ près~? Comparer la valeur de $f(7)$ et la valeur absolue du plus grand terme de la série, quelle est la perte de précision relative si on effectue les calculs en virgule flottante~? Combien de chiffres significatifs faut-il utiliser pour assurer une précision finale de 16 chiffres en base 10~? Calculer le développement asymptotique en l'infini et déterminer un encadrement de $f(7)$ par ce développement. Combien de termes faut-il calculer pour déterminer $f(10)$ à $1e-16$ près par le développement asymptotique et par le développement en séries~? Quelle est la meilleure méthode pour calculer $f(10)$~? \subsubsection{Recherche de solutions d'\'equations diff\'erentielles} On peut aussi appliquer les techniques ci-dessus pour calculer des solutions de certaines \'equations diff\'erentielles dont les solutions ne s'expriment pas \`a l'aide des fonctions usuelles, on remplace dans l'\'equation la fonction inconnue par son d\'eveloppement en s\'eries et on cherche une relation de r\'ecurrence entre $a_{n+1}$ et $a_n$. Si on arrive \`a montrer par exemple qu'il y a une solution ayant un d\'eveloppement altern\'ee, ou plus g\'en\'erallement, si on a une majoration $|a_{n+1}/a_n| |x|$, \[ |a_{2k} x^{2k}| < |a_{2k-2} x^{2k-2}| \] donc la s\'erie $\sum_{k=0}^\infty a_k x^k$ est altern\'ee \`a partir du rang partie enti\`ere de $|x|$ plus un. Donc elle converge pour tout $x$ (le rayon de convergence de $y$ est $+\infty$) et le reste de la somme jusqu'\`a l'ordre $2n$ est inf\'erieur en valeur absolue \`a~: \[ |R_{2n}(x)| \leq |a_{2n+2} x^{2n+2}| \] Par exemple, pour avoir une valeur approch\'ee \`a $1e-10$ pr\`es de $y(x)$ pour $a_0=1$ et $|x|\leq 1$, on calcule $y=\sum_{k=0}^{2n} a_k x^k $, on s'arr\^ete au rang $n$ tel que \[ |a_{2n+2} x^{2n+2}| \leq |a_{2n+2}| \leq 10^{-10} \] On remarque que~: \[ a_{2n} = \frac{(-1)^n}{2^2 4^2 ... (2n)^2} = \frac{(-1)^n}{2^{2n} n!^2} \] donc $n=7$ convient. Pour $m \neq 0$, on peut faire un raisonnement analogue (les calculs sont un peu plus compliqu\'es). On a ainsi trouv\'e une solution $y_0$ de l'\'equation diff\'erentielle de d\'epart dont on peut facilement calculer une valeur approch\'ee (aussi facilement que par exemple la fonction sinus pour $|x| \leq 1$), on peut alors trouver toutes les solutions de l'\'equation diff\'erentielle (en posant $y=y_0 z$ et en cherchant $z$). {\bf Exercice~:} faire de m\^eme pour les solutions de $y'{'}-xy=0$ (fonctions de Airy). \subsection{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)=\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^{-t}$ (qui admet $-e^{-t}$ comme primitive, cette primitive ayant une limite en $+\infty$). Pour trouver le développement asymptotique de $f$ 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 \begin{eqnarray*} 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{eqnarray*} où \begin{equation} \label{eq:Eiinf} 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 \end{equation} Le développement en séries est divergent puisque pour $x>0$ fixé et $n$ tendant vers l'infini \[ \lim_{n\rightarrow +\infty} \frac{n!}{x^{n+1}} = +\infty\] mais si $x$ est grand, au début la série semble converger, de manière très rapide~: \[ \frac{1}{x} >> \frac{1}{x^2} >> \frac{2}{x^3} \] On peut utiliser $S(x)$ comme valeur approchée de $f(x)$ pour $x$ grand si on sait majorer $R(x)$ par un nombre suffisamment petit. On a \[ | R(x) | \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\'e par $e^{-x}/x$. Pour $x$ fixé assez grand, il faut donc trouver un rang $n$, s'il en existe un, tel que $(n+1)!/x^{n+1}<\epsilon$ où $\epsilon$ est la précision relative que l'on s'est fixée. Par exemple, si $x\geq 100$, $n=11$ convient pour $\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\'ecision et peu de calculs, mais contrairement aux s\'eries enti\`eres, il n'est pas possible d'am\'eliorer cette pr\'ecision de mani\`ere arbitraire en poussant le d\'eveloppement plus loin, il y a une pr\'ecision maximale possible (qui d\'epend de $x$). Ce type de d\'eveloppement asymptotique peut \^etre effectu\'e pour d'autres fonctions du m\^eme type, par exemple \[ \int_x^{+\infty} e^{-t^2} \ dt, \quad \int_x^{+\infty} \frac{\sin(t)}{t} \ dt, \quad ... \] {\bf Digression: calcul approché de la constante d'Euler $\gamma$}\\ On peut montrer que \begin{equation} \label{eq:def_gamma} \lim_{n\rightarrow +\infty} u_n, \quad u_n=\sum_{k=1}^{n}\frac{1}{k} - \ln(n) \end{equation} existe (par exemple en cherchant un \'equivalent de $u_{n+1}-u_n$ qui vaut $\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 $x\rightarrow 0$, $f(x)$ admet une singularité en $-\ln(x)$, plus précisément $f(x)+\ln(x)$ admet un développement en séries (de rayon de convergence $+\infty$), car~: \begin{eqnarray*} 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{eqnarray*} Que vaut la constante du membre de droite~: \[ 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=-\gamma$ (voir plus bas une démonstration condensée) et donc~: \begin{equation} \label{eq:gamma} \gamma= \int_0^{x} \frac{1-e^{-t}}{t} \ dt -f(x)-\ln(x) \end{equation} Pour obtenir une valeur approchée de $\gamma$, il suffit donc de prendre un $x$ assez grand pour pouvoir calculer $f(x)$ par son développement asymptotique à la précision requise, puis de calculer l'intégrale du membre de droite par le développement en séries en $x=0$ (en utilisant une précision intermédiaire plus grande puisque ce développement en séries va sembler diverger au début avant de converger pour $n$ suffisamment grand). Par exemple, on pose $x=13$, on calcule $f(13)$ par (\ref{eq:Eiinf}) avec $n=13$ (qui correspond au moment o\`u le terme g\'en\'eral de la s\'erie est minimum puisque le rapport de deux termes successifs est en $n/x$) et une erreur absolue inf\'erieure \`a $e^{-13} 13!/13^{14}=4e-12$ \begin{center} $f(13) \approx$ \verb|exp(-13)*sum((-1)^n*n!/13.^(n+1),n=0..13)| \end{center} puis on remplace dans (\ref{eq:gamma}), avec \[ \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\'ee, en faisant la somme jusqu'au rang 49 (pour lequel le terme g\'en\'eral est de l'ordre de \verb|1e-12|), le reste de cette somme $R_{50}$ est positif et est inf\'erieur \`a {\tt (-1)\verb|^|50*13.\verb|^|51/51/51!)} qui est de l'ordre de \verb|8e-12| \begin{center} \verb|evalf(sum((-1)^n*13^(n+1)/(n+1)/(n+1)!,n=0..49))| \end{center} La somme argument de \verb|evalf| \'etant exacte, il n'y a pas de probl\`emes de perte de pr\'ecision, on peut aussi faire les calculs interm\'ediaires en arithm\'etique approch\'ee, on doit alors prendre 4 chiffres significatifs de plus pour tenir compte de la valeur du plus grand terme somm\'e dans la s\'erie, terme que l'on d\'etermine par exemple par \begin{center} {\tt seq(13.\verb|^|(n+1)/(n+1)/(n+1)!,n=0..20)} \end{center} ce terme vaut \verb|13^11/11/11!| soit 4000 environ) \begin{center} \verb|Digits:=16; sum((-1)^n*13.^(n+1)/(n+1)/(n+1)!,n=0..49)| \end{center} On obtient finalement comme valeur approch\'ee de $\gamma$ \begin{center} \verb|-exp(-13)*sum((-1)^n*n!/13.^(n+1),n=0..13)-ln(13)+|\\ \verb| sum((-1)^n*13^(n+1)/(n+1)/(n+1)!,n=0..49)| \end{center} soit \verb|0.577215664897| avec une erreur inf\'erieure \`a \verb|1.2e-11|. Bien entendu, cette m\'ethode est surtout int\'eressante si on veut calculer un grand nombre de d\'ecimales de la constante d'Euler, sinon on peut par exemple appliquer la m\'ethode d'acc\'el\'eration de Richardson \`a la suite convergente (\ref{eq:def_gamma}) qui d\'efinit $\gamma$ ou d'autres m\'ethodes d'acc\'el\'eration (en transformant par exemple la s\'erie en s\'erie altern\'ee). On calcule alors de deux mani\`eres diff\'erentes $f(x)$ pour $x$ plus grand (d\'etermin\'e par la pr\'ecision qu'on peut obtenir par le d\'eveloppement aymptotique de $f$). 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 $f$) et l'égalité (dont un sch\'ema de preuve est aussi donn\'e plus bas) \begin{equation} \label{eq:Siinf} \int_0^{+\infty} \frac{\sin(t)}{t} \ dt = \frac{\pi}{2} \end{equation} {\bf Calcul de $C$ (et preuve de (\ref{eq:Siinf}))}:\\ Pour cela on effectue une intégration par parties, cette fois en intégrant $1/t$ et en dérivant l'exponentielle (moins 1 dans la première intégrale). \begin{eqnarray*} 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{eqnarray*} Pour calculer cette intégrale, on utilise l'égalité (qui se démontre par récurrence en faisant une intégration par parties)~: \[ 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 $t$ et on dérive l'autre terme, puis $t^2/2$, etc. \begin{eqnarray*} 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{eqnarray*} où \[ I_n=\int_0^{+\infty} \frac{t^n}{n!} e^{-t} (\ln(t)-\ln(n)) \ dt \] Pour déterminer $I_n$ on fait le changement de variables $t=nu$ \begin{eqnarray*} 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{eqnarray*} Or en faisant le même changement de variables $t=nu$~: \[ n!= \int_0^{+\infty}t^n e^{-t} \ dt = n^{n+1} \int_0^{+\infty} e^{n(ln(u)-u)} \ du \] Donc \[ I_n= \frac{\int_0^{+\infty} e^{n(ln(u)-u)} \ln(u) \ du} {\int_0^{+\infty} e^{n(ln(u)-u)} \ du} \] Lorsque $n$ tend vers l'infini, on peut montrer que $I_n \rightarrow 0$, en effet les intégrales sont équivalentes à leur valeur sur un petit intervalle autour de $u=1$, point où l'argument de l'exponentielle est maximal, et comme l'intégrand du numérateur a une amplitude $\ln(u)$ qui s'annule en $u=1$, il devient négligeable devant le dénominateur. Finalement on a bien $C=-\gamma$. On peut remarquer qu'en faisant le m\^eme calcul que $C$ mais en remplacant $e^{-t}$ par $e^{-\alpha t}$ pour $\Re(\alpha)>0$, donne $\lim I_n=-\ln(\alpha)$ (car le point critique o\`u la d\'eriv\'ee de la phase s'annule est alors $1/\alpha$). Ceci peut aussi se v\'erifier pour $\alpha$ r\'eel en faisant le changement de variables $\alpha t=u$ \[ \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$, $-\ln(\alpha)$ tend vers $\ln(i)=i\frac{\pi}{2}$ et on obtient \[ \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 (\ref{eq:Siinf}), et la partie r\'eelle une autre identit\'e sur $\gamma$ faisant intervenir la fonction cosinus int\'egral. \section{Polynômes~: arithm\'etique, factorisation, interpolation} \subsection{Arithm\'etique des polynomes: B\'ezout et applications} On consid\`ere les polyn\^omes \`a une variable \`a coefficients dans $\R$ ou $\C$ ou $\Q$. Les algorithmes de base d\'ej\`a \'evoqu\'es sont l'\'evaluation en un point (m\'ethode de Horner), l'addition, la soustraction, la multiplication et la division euclidienne de $A$ par $B \neq 0$~: \[ A = B Q + R, \quad \mbox{deg}(R)< \mbox{deg}(B) \] A l'aide de la division euclidienne, on peut calculer le PGCD de deux polyn\^omes par l'algorithme d'Euclide. Nous allons pr\'esenter l'algorithme d'Euclide \'etendu (ou de B\'ezout) \begin{theorem} \'Etant donn\'es 2 polyn\^omes $A$ et $B$, il existe deux polyn\^omes $U$, $V$ tels que \[ AU + B V = \mbox{pgcd} (A,B) , \quad \mbox{deg}(U) < \mbox{deg}(B), \mbox{deg}(V) < \mbox{deg}(A) \] \end{theorem} Algorithme~: \index{Bezout}\\ On construit en fait 3 suites $(U_n)$, $(V_n)$ et $(R_n)$ telles que~: \[ AU_n + B V_n = R_n \] \begin{itemize} \item on initialise $U_0=1, V_0=0, R_0=A$ et $U_1=0, V_1=1, R_1=B$ \item on calcule les indices $n+2$ en fonction de $n$ et $n+1$ en effectuant la division euclidienne de $R_n$ par $R_{n+1}$ \[ R_n = Q_n R_{n+1} + R_{n+2}, \quad U_{n+2}= U_n - Q_n U_{n+1}, V_{n+2} = V_n - Q_n V_{n+1} \] \item on s'arr\^ete au dernier reste non nul \end{itemize} Exemple~:\\ $A=x^3-1, B=x^2+1$, les rangs 0 et 1 sont donnés ci-dessus. Au rang 2, $Q_0$ est le quotient euclidien de $A$ par $B$ (fonction \verb|quo|) donc $x$, d'où \[ U_2=1, V_2= -x, R_2=-x-1\] Puis on divise $x^2+1$ par $-x-1$, quotient $-x+1$, donc \[ 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\'ecurrence que la relation $AU_n + B V_n = R_n$ est conserv\'ee. Comme $R_n$ est la suite des restes, le dernier reste non nul est bien le pgcd de $A$ et $B$. D'autre part, examinons les degr\'es des $V_k$. Supposons que deg$(A) \geq $ deg$(B)$ (sinon on \'echange $A$ et $B$). Au rang $n=0$, $V_0=0$ donc $V_2=-Q_0 V_1$, aux rangs suivants le degr\'e de $Q_n$ est non nul (car le degr\'e de $R_{n+1}$ est strictement inf\'erieur au degr\'e de $R_n$) On montre donc par r\'ecurrence que la suite des degr\'es de $V_n$ est croissante et que~: \[ \mbox{deg} (V_{n+2}) = \mbox{deg}(Q_n)+\mbox{deg}(V_{n+1}) \] Comme deg$(Q_n)$=deg$(R_n)$-deg$(R_{n+1})$, on en d\'eduit que \[ \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+2$ est le rang du dernier reste non nul, $V_{n+2}=V$ et deg$V$=deg$A$-deg$R_{n+1}$ est donc strictement inf\'erieur au degr\'e de $A$ (car $R_{n+1}$, l'avant-dernier reste non nul, est de degr\'e plus grand ou \'egal \`a 1). On en d\'eduit enfin que le degr\'e de $U$ est strictement inf\'erieur au degr\'e de $B$, car $AU=R-BV$, le degr\'e de $BV$ est strictement inf\'erieur \`a celui de $B$ plus celui de $A$. L'identit\'e de B\'ezout permet de r\'esoudre plus g\'en\'erallement une \'equation du type \[ Au+Bv=C\] o\`u $A,B,C$ sont trois polyn\^omes donn\'es, \`a condition que $C$ soit divisible par le pgcd de $A$ et $B$. L'ensemble des solutions s'obtient \`a partir d'une solution particuli\`ere $U,V$ de B\'ezout, notons $c=C/\mbox{gcd}(A,B)$, on a alors \[ A(cU) + B(cV)=c\ \mbox{gcd}(A,B) = C \] et l'ensemble des solutions est donn\'e par $u=cU-PB, v=cV+PA$ o\`u $P$ est un polyn\^ome quelconque. Si le degr\'e de $C$ est plus petit que le degr\'e de $A$ plus le degr\'e de $B$, il existe une solution ``priviligi\'ee'', on prend pour $u$ le reste de la division euclidienne de $cU$ par $B$, $v$ est alors le reste de la division euclidienne de $cV$ par $A$ pour des raisons de degr\'e. Exemple~: si on veut résoudre \[ (x^3-1)u+(x^2+1)v=2x^2\] on multiplie $U=x-1$ et $V=1+x-x^2$ par $x^2$ ce qui donne une solution \[ 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), \quad v-P(x^3-1) \] et la solution priviligiée (de degrés minimaux) est \[ -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\'e de B\'ezout intervient dans de nombreux probl\`emes en particulier la d\'ecomposition en \'el\'ements simples d'une fraction rationnelle. Si le d\'enominateur $D$ d'une fraction se factorise en produit de 2 facteurs $D=AB$ premiers entre eux, alors il existe deux polyn\^omes $u$ et $v$ tels que $N=Au+Bv$, donc \[ \frac{N}{D} = \frac{Au+Bv} {AB} = \frac{u}{B} + \frac{v}{A} \] Si de plus $N/D$ est une fraction propre (degr\'e de $N$ plus petit que celui de $D$), alors $u/B$ et $v/A$ sont encore des fractions propres (en calculant le reste de la division euclidienne pour $u$ et $v$ comme expliqu\'e ci-dessus). Par exemple~: \[ \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 \begin{itemize} \item le calcul de primitive de fraction rationnelles (et tout ce qui s'y ram\`ene), par exemple \[ \int \frac{2x^2}{(x^3-1)(x^2+1)} = = \int \frac{-x+1}{x^2+1} + \int \frac{x^2-x+1}{x^3-1} \] Puis on fait apparaitre la dérivée du dénominateur au numérateur pour éliminer les $x$, $2x=(x^2+1)'$ \begin{eqnarray*} \int \frac{-x+1}{x^2+1} &= &-\frac{1}{2} \int \frac{(x^2+1)'}{x^2+1} + \int \frac{1}{x^2+1} + \int \frac{x^2-x+1}{x^3-1} \\ &=& -\frac{1}{2} \ln(x^2+1) + \arctan(x)+\int \frac{x^2-x+1}{x^3-1} \end{eqnarray*} pour faire le calcul complet, il faut aussi décomposer la fraction restante (exercice!) \item Le calcul de la fonction exponentielle (\`a nouveau). Au lieu d'utiliser $T$ le d\'eveloppement de Taylor en 0 par exemple \`a l'ordre 10, on cherche une fraction rationnelle $N/D$ ayant le m\^eme d\'eveloppement de Taylor que l'exponentielle en 0 avec degr\'e de $N$ et de $D$ major\'es par 5. Pour trouver $N$ et $D$ on multiplie la condition $N/D=T+O(x^{11})$ par $D$ ce qui donne \[ N=DT+O(x^{11})=DT+Px^{11}\] on applique l'algorithme de B\'ezout aux polyn\^omes $x^{11}$ et $T$ en s'arr\^etant pr\'ematur\'ement, lorsque le reste est de degr\'e 5, on montre alors que le reste est $N$ et le coefficient de B\'ezout de $T$ est $D$. On peut alors montrer que l'approximation est un peu meilleure, et n\'ecessite moins d'op\'erations (il y a une certaine sym\'etrie entre les termes de $N$ et $D$). \item le calcul de transform\'ee de Laplace inverse de fractions rationnelles, l'idée est la même, sauf qu'on remplace l'intégrale par la transformée de Laplace inverse (et les formules donnant la transformée inverse de $1/(x-p)$, $1/(x^2+p^2)$, $p/(x^2+p^2)$ respectivement $\exp(px), \sin(xp)/p, \cos(px)$) (calcul non exigible à l'examen) \item le calcul du terme d'ordre $n$ du d\'eveloppement de Taylor en 0 d'une fraction rationnelle. On décompose, et on se ramène à des séries dont le terme général est connu, comme $(a+x)^{-n}$. Par exemple pour connaitre le développement de $1/(x^2-3x+2)$, on factorise le dénominateur $1/((x-1)(x-2))$, on décompose \[ \frac{1}{(x-1)(x-2)} = \frac{-1}{x-1} + \frac{1}{x-2} = \frac{1}{1-x} - \frac{1}{2} \frac{1}{1-\frac{x}{2}}\] et on développe, le terme d'ordre $n$ est donc $1-(1/2)^ {n+1}$. %\item Le calcul de sommes de séries de certaines fractions rationnelles. \end{itemize} Il faut n\'eanmoins savoir factoriser un polyn\^ome, ce dont nous parlerons dans la section suivante. Exercice~: Calculer l'intégrale \[ \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 $n$ 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'\'elimination dans les syst\`emes polynomiaux, par exemple consid\'erons le syst\`eme de 2 \'equations \`a 2 inconnues (intersection d'une ellipse et d'un cercle)~: \[ x^2+y^2-9=0, x^2+2y^2-2xy-7=0\] En calculant les coefficients de B\'ezout des 2 polyn\^omes en $x$ $x^2+y^2-9$ et $x^2+2y^2-2xy-7$ et en multipliant au besoin par le PPCM (plus grand commun multiple) des d\'enominateurs, on obtient \`a droite de l'\'equation de B\'ezout un polyn\^ome ne d\'ependant que de $y$ et qui s'annule aussi aux solutions du syst\`eme, on peut alors r\'esoudre en $y$ (en factorisant) puis en $x$. Ici par exemple ce polynome est $5y^4-32y^2+4$. Cette m\'ethode se syst\'ematise, le polynome obtenu par \'elimination d'une variable est appel\'e r\'esultant. %Liens: %\verb|formation.etud.u-psud.fr/maths/mathsenligne.pdf|\\ %\verb|http://wims.unice.fr/wims/fr_tool~arithmetic~bezout.html| \subsection{Factorisation des polyn\^omes} Soit $P$ un polyn\^ome de degr\'e non nul. Factoriser $P$ n'a pas une signification unique, tout d\'epend d'une part si on veut une factorisation exacte ou approch\'ee, et d'autre part quels seront les types des coefficients de la factorisation (complexes, r\'eels, entiers). \subsubsection{Multiplicit\'e des racines.} \index{factorisation} \index{multiplicite} \index{racine} On dit que $r$ est une racine de multiplicit\'e $k$ de $P$ si $P(x)=(x-r)^k Q$ et $Q(r)\neq 0$. En faisant le d\'eveloppement de Taylor de $P$ en $r$ \`a l'ordre degr\'e de $P$, on voit que cela \'equivaut \`a~: \[ P(r)=P'(r)=...=P^{[k-1]}(r)=0, \quad P^{[k]}(r) \neq 0 \] En particulier si $P(r)=0$, on peut factoriser $P$ par $X-r$. On peut donc d\'etecter les racines de multiplicit\'e sup\'erieure \`a 1 en cherchant un facteur commun \`a $P$ et $P'$, en effet $x-r$ divisera $P$ et $P'$. \begin{theorem} Si $P$ et $P'$ sont premiers entre eux (pgcd = 1), alors les racines de $P$ sont simples (de multiplicit\'e 1). \end{theorem} Il existe un algorithme (d\^u \`a Yun) qui permet d'\'ecrire un polynome quelconque comme produit de polyn\^omes dont les racines sont simples en effectuant uniquement des calculs de PGCD de polynomes. \index{squarefree} \begin{verbatim} 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)); }; } \end{verbatim} L'instruction \verb|squarefree| ou \'equivalente de votre logiciel de calcul formel effectue cette d\'ecomposition. \subsubsection{Factorisation dans $\C$.} Reste maintenant \`a trouver des racines! On a le~: \begin{theorem} (d'Alembert) \index{factorisation}\\ Soit $P$ un polynome de degr\'e non nul, alors $P$ admet au moins une racine complexe. \end{theorem} On peut alors factoriser $P$ par $X-r$ si $r$ est la racine, et recommencer avec le quotient, d'o\`u le corollaire. \begin{theorem} Soit $P$ un polynome de degr\'e $n$ non nul, alors $P$ admet $n$ racines complexes (compt\'ees avec multiplicit\'e) $x_1,...,x_n$, on a donc~: \begin{equation} \label{eq:factorisation} P(X)=a_n \prod_{j=1}^n (X-x_j) \end{equation} o\`u $a_n$ est le coefficient dominant de $P$. \end{theorem} D\'emonstration du th\'eor\`eme de d'Alembert~:\\ On va montrer que le minimum de la valeur absolue de $P$ est atteint en un nombre complexe puis que ce minimum est forc\'ement nul. Soit \[ P(x)=a_n x^n + ... + a_0, \quad a_n \neq 0 \] Lorsque $|x|$ tend vers l'infini, $|P(x)|$ tend vers l'infini, en effet \[ 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\'ecis\'ement il existe $R>0$ tel que si $|x|>R$ alors $|P(x)|>|a_n| |x|^n/2$. Quitte \`a augmenter $R$ on peut donc supposer que $|P(x)|>|P(0)|$ si $|x|>R$, donc il existe un complexe $x_0$ qui r\'ealise le minimum de $|P|$ sur $\C$ (ce minimum est en fait le minimum pour $|x|\leq R$). On va montrer par l'absurde que ce minimum est nul (donc que $x_0$ est la racine cherch\'ee). Supposons donc que $P(x_0) \neq 0$. On fait le d\'eveloppement de Taylor de $P$ en $x_0$ \`a l'ordre $n$=degr\'e de $P$, donc le d\'eveloppement n'a pas de reste~: \[ P(x) - P(x_0) = (x-x_0) P'(x_0) + .. + (x-x_0)^n \frac{P^{[n]}(x_0)}{n!}\] Comme $P$ n'est pas constant, l'un des termes du membre de droite est non nul, soit $k$ l'indice du premier terme non nul, on a alors~: \[ P(x) = P(x_0) + (x-x_0)^k \frac{P^{[k]}(x_0)}{k!} + o((x-x_0)^k) \] Comme $P(x_0) \neq 0$, on peut le factoriser en~: \[ 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+t w$ o\`u $w$ est une racine $k$-i\`eme (cela existe dans $\C$) de \[ \left(\frac{-P^{[k]}(x_0)}{P(x_0) k!} \right)^{-1} \] on a alors~: \[ P(x)=P(x_0)(1-t^k+o(t^{k+1}))\] lorsque $t$ est positif, suffisamment petit, on a $0 < 1-t^k+o(t^{k+1} < 1$, donc $|P(x)|<|P(x_0)|$, ce qui est absurde ($x_0$ r\'ealisant le minimum de $P$ sur $\C$). {\bf Remarque}~:\\ Si on d\'eveloppe la relation (\ref{eq:factorisation}), on obtient des relations entre les coefficients du polynome et les racines, par exemple~: \[ a_{n-1}=a_n \sum{j=1}^n (-x_j), ...,\quad a_0 = a_n \prod_{j=1}^n (-x_j), \] \subsubsection{Calcul approch\'e des racines complexes simples} \index{factorisation} \index{racine} La section pr\'ec\'edente nous a montr\'e qu'on pouvait se ramener \`a la recherche de racines simples, ce qui donne envie d'essayer la m\'ethode de Newton. On a malheureusement rarement la possibilit\'e de pouvoir d\'emontrer qu'\`a partir d'une valeur initiale donn\'ee, la m\'ethode de Newton converge, parce que les racines peuvent \^etre complexes, et m\^eme si elles sont r\'eelles, on n'a pas forc\'ement de r\'esultat sur la convexit\'e du polyn\^ome (cf. cependant une application des suites de Sturm dans la section suivante qui permet de connaitre le signe de $P'{'}$ sur un intervalle sans le factoriser). On effectue donc souvent des it\'erations de Newton, en partant de 0.0, en esp\'erant s'approcher suffisamment d'une racine pour que le th\'eor\`eme de convergence th\'eorique s'applique. On se fixe un nombre maximal d'it\'erations, si on le d\'epasse on prend alors une valeur initiale al\'eatoire complexe et on recommence. Une fois une racine d\'etermin\'ee, on l'\'elimine en calculant le quotient euclidien $Q$ de $P$ par $X-r$ (par l'algorithme de Horner), puis on calcule les racines du quotient $Q$ (qui sont des racines de $P$). Un probl\`eme pratique apparait alors, c'est que $r$ n'est pas exact donc le quotient $Q$ non plus, au fur et \`a mesure du calcul des racines de $P$, on perd de plus en plus de pr\'ecision. Il existe une am\'elioration simple, si $r'$ est une racine approch\'ee de $Q$, alors elle est racine approch\'ee de $P$ et on a toutes les chances qu'elle soit suffisamment proche d'une racine de $P$ pour que le th\'eor\`eme s'applique, on effectue alors 1 ou 2 it\'erations de Newton avec $r'$ mais pour $P$ (et non $Q$) afin d'am\'eliorer sa pr\'ecision comme racine de $P$. \subsubsection{Factorisation dans $\R$, localisation des racines} \index{factorisation} Pour factoriser un polyn\^ome \`a coefficients r\'eels, on commence par le factoriser dans $\C$. On observe ensuite que si $r$ est une racine complexe non r\'eelle de $P$, alors son conjugu\'e l'est aussi (il suffit de prendre le conjugue de la relation $P(r)=0$) et avec la m\^eme multiplicit\'e (les d\'eriv\'ees successives de $P$ \'etant aussi \`a coefficients r\'eels). On regroupe alors les facteurs correspondant \`a des racines complexes conjugu\'ees~: \[ (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~: \begin{theorem} La factorisation d'un polyn\^ome \`a coefficients r\'eels sur $\R$ donne un produit de facteurs de degr\'e 1 (correspondant \`a des racines r\'eelles) et de degr\'e 2 (correspondant \`a des paires de racines complexes conjugu\'ees) \end{theorem} Il existe un algorithme utilisant l'algorithme de calcul du PGCD de $P$ et $P'$ qui permet de d\'eterminer le nombre de racines r\'eelles d'un polyn\^ome $P$ sans racine multiple sur $\R$ ou dans un intervalle de $R$. \begin{theorem} \index{Sturm} On d\'efinit la suite de polyn\^omes $A_0=P, A_1=P', ..., A_k,0$ en prenant l'oppos\'e du reste de la division euclidienne des deux pr\'ec\'dents~: \begin{equation} \label{eq:sturm} A_{i} = A_{i+1} Q_{i+2} - A_{i+2} \end{equation} Soit $A_k$, le dernier reste non nul, c'est un polyn\^ome constant puisque $P$ n'a pas de racine multiple. On d\'efinit $s(a)$ comme \'etant le nombre de changements de signes de la suite $A_i(a)$ en ignorant les 0. Alors le nombre de racines r\'eelles de $A_0=P$ sur l'intervalle $]a,b]$ est \'egal \`a $s(a)-s(b)$. \end{theorem} {\bf Exemple}~:\\ Quel est le nombre de racines r\'eelles de $P=x^3+x+1$ sur $[-2,2]$? sur $[0,2]$?\\ On a donc \[ 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=-2$ on obtient la suite $-9,13,1/3,-31/4$ (2 changements de signe), en $x=2$ on obtient la suite $11,13,-7/3,-31/4$ (1 changement de signe), il y a donc 1 racine r\'eelle entre -2 et 2. En $x=0$ on obtient la suite $1,1,-1,-31/4$ (1 changement de signe) donc la racine r\'eelle est entre -2 et 0. {\bf Preuve}\\ On consid\'ere 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 (\ref{eq:sturm}), or $A_k$ est constant non nul). Elle ne peut pas non plus contenir ...,+,0,+,... ni ...,-,0,-,... \`a cause de la convention de signe sur les restes de (\ref{eq:sturm}). Donc si $b$ est une racine de $A_i$ pour $00$ donc il existe une racine $r>1$ telle que $P'(r)>0$, toute valeur de départ de Newton supérieure à $r$ assure la convergence. {\bf Remarque}~:\\ On peut aussi d\'eterminer les racines r\'eelles d'un polyn\^ome \`a coefficients rationnels en faisant uniquement des calculs exacts par dichotomie. Cette m\'ethode de localisation des racines r\'eelles se g\'en\'eralise d'ailleurs au cas complexe. On peut ainsi d\'eterminer les racines complexes d'un polyn\^ome \`a coefficients complexes rationnels de mani\`ere d\'eterministe \`a la pr\'ecision voulue (cf. Eisermann). \subsubsection{Factorisation exacte} \index{factorisation} Soit $P$ un polyn\^ome \`a coefficients entiers. Lorsqu'on demande \`a un logiciel de calcul formel de factoriser $P$, par d\'efaut il ne calcule pas les racines complexes approch\'ees, mais renvoie une factorisation {\bf exacte}, sous forme de produit de facteurs \`a coefficients {\bf entiers}. Les degr\'es des facteurs peuvent \^etre plus grand que 2. Par exemple $x^4+x+1$ ne peut pas \^etre factoris\'e en produit de polyn\^omes \`a coefficients entiers (bien qu'il ait 2 facteurs de degr\'e 2 dans $\R$ et 4 de degr\'e 1 dans $\C$). \index{racines rationnelles} Commencons par une m\'ethode simple de calcul des racines rationnelles de $P$ (les racines rationnelles correspondent \`a des facteurs entiers de degr\'e 1 de la forme $qX-p$ de $P$). Soit $x=p/q$ une racine rationnelle \'ecrite sous forme de fraction irr\'eductible de $P=a_n X^n+...+a_0$, on a alors \[ 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_n p^{n-1} + a_{n-1} p^{n-1}q+...+a_1 q^{n-1}) = - a_0 q^n \] et $p$ divise donc $a_0 q^n$. Comme $p/q$ est irr\'eductible, cela entraine que $p$ divise $a_0$. De m\^eme $q$ divise $a_n$. Il suffit donc de tester quelles sont les racines de $P$ parmi toutes les fractions irr\'eductibles de la forme un diviseur de $a_0$ sur un diviseur de $a_n$ (attention \`a ne pas oublier les diviseurs n\'egatifs!). Exemple: racines rationnelles de $2x^2+3x+1=0$. On a $p$ divise $1$ donc vaut 1 ou -1, $q$ divise $2$ donc vaut 1 ou 2. On teste donc 1, -1, 1/2, -1/2. On obtient ici la factorisation compl\`ete du polynome (les racines sont -1 et -1/2) \[ 2x^2+3x+1=2(x+1)(x+1/2) \] {\bf Remarques}~: \begin{itemize} \item Pour un polynome al\'eatoire, on ne trouvera aucune racine rationnelle. \item Cette m\'ethode n'est pas tr\`es efficace, car factoriser un entier peut \^etre long, le nombre de tests peut \^etre tr\`es grand (si $a_n$ et $a_0$ ont beaucoup de facteurs), les logiciels de calcul formel utilisent des m\'ethodes appel\'ees $p$-adiques pour trouver les racines rationnelles d'un polynome (on calcule d'abord les racines de $P$ modulo $p$ puis modulo $p^k$ pour $k$ assez grand). On pourrait aussi penser à calculer les racines complexes approchées et voir si en multipliant par $a_n$ on est proche d'un entier, on testerait alors le rationnel correspondant. \end{itemize} Pour d\'eterminer les facteurs \`a coefficients entiers de plus grand degr\'e, il n'existe pas de m\'ethode aussi simple. On peut calculer des valeurs approch\'ees des racines complexes et essayer de cr\'eer des paquets de racines complexes, puis tester si $a_n\prod_{r \in \mbox{paquet}} (X-r)$ est \`a coefficients entier (aux erreurs d'arrondi pr\`es). Par exemple si on calcule les racines complexes approch\'ees de $x^6+2x^3-x^2+1$, on pourra composer un facteur de degr\'e 3 \`a coefficients entiers en rassemblant les racines de $x^3+x+1$. Les logiciels de calcul formel utilisent des algorithmes modulaires et $p$-adiques (consistant \`a factoriser le polynome modulo $p$). \subsection{Approximation polynomiale} \'Etant donn\'e la facilit\'e de manipulation qu'apportent les polynomes, on peut chercher \`a approcher une fonction par un polyn\^ome. La m\'ethode la plus naturelle consiste \`a chercher un polyn\^ome de degr\'e le plus petit possible \'egal \`a la fonction en certains points $x_0,...,x_n$ et \`a trouver une majoration de la diff\'erence entre la fonction et le polyn\^ome. Le polynome interpolateur de Lagrange r\'epond \`a cette question. Soit donc $x_0,...,x_n$ des r\'eels distincts et $y_0,...,y_n$ les valeurs de la fonction \`a approcher en ces points (on posera $y_j=f(x_j)$ pour approcher la fonction $f$). On cherche donc $P$ tel que $P(x_j)=y_i$ pour $j \in [0,n]$. Commencons par voir s'il y a beaucoup de solutions. Soit $P$ et $Q$ deux solutions distinctes du probl\`eme, alors $P-Q$ est non nul et va s'annuler en $x_0, ...,x_n$ donc poss\`ede $n+1$ racines donc est de degr\'e $n+1$ au moins. R\'eciproquement, si on ajoute \`a $P$ un multiple du polynome $A=\prod_{j=0}^n (X-x_j)$, on obtient une autre solution. Toutes les solutions se d\'eduisent donc d'une solution particuli\`ere en y ajoutant un polynome de degr\'e au moins $n+1$ multiple de $A$. Nous allons maintenant construire une solution particuli\`ere de degr\'e au plus $n$. Si $n=0$, on prend $P=x_0$ constant. On proc\`ede ensuite par r\'ecurrence. Pour construire le polyn\^ome correspondant \`a $x_0,...,x_{n+1}$ on part du polyno\^ome $P_n$ correspondant \`a $x_0,...,x_{n}$ et on lui ajoute un multiple r\'eel de $A$ \[ 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_j$ pour $j=0,..n$, on calcule maintenant $ \alpha_{n+1}$ pour que $P_{n+1}(x_{n+1})=y_{n+1}$. En remplacant avec l'expression de $P_{n+1}$ ci-dessus, on obtient \[ P_n(x_{n+1})+ \alpha_{n+1} \prod_{j=0}^n (x_{n+1}-x_j) = y_{n+1} \] Comme tous les $x_j$ sont distincts, il existe une solution unique~: \[ \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\'e le~: \begin{theorem} \index{lagrange} \index{interpolation} Soit $n+1$ r\'eels distincts $x_0,...,x_n$ et $n+1$ r\'eels quelconques $y_0,...,y_n$. Il existe un unique polyn\^ome $P$ de degr\'e inf\'erieur ou \'egal \`a $n$, appel\'e polynome de Lagrange, tel que~: \[ P(x_i)=y_i\] \end{theorem} Exemple~: d\'eterminons le polynome de degr\'e inf\'erieur ou \'egal \`a 2 tel que $P(0)=1, P(1)=2, P(2)=1$. On commence par $P_0=1$. Puis on pose $P_1=P_0+ \alpha_{1}X=1+ \alpha_{1}X$. Comme $P(1)=2=1+ \alpha_{1}$ on en tire $ \alpha_{1}=1$ donc $P_1=1+X$. Puis on pose $P_2=P_1+ \alpha_{2}X(X-1)$, on a $P_2(2)=3+2 \alpha_{2}=1$ donc $ \alpha_{2}=-1$, finalement $P_2=1+X-X(X-1)$. Reste \`a estimer l'\'ecart entre une fonction et son polynome interpolateur, on a le~: \begin{theorem} \index{lagrange} Soit $f$ une fonction $n+1$ fois d\'erivable sur un intervalle $I=[a,b]$ de $\R$, $x_0,...,x_n$ des r\'eels distincts de $I$. Soit $P$ le polynome de Lagrange donn\'e par les $x_j$ et $y_j=f(x_j)$. Pour tout r\'eel $x \in I$, il existe un r\'eel $\xi_x \in [a,b]$ (qui d\'epend de $x$) tel que~: \begin{equation} \label{eq:lagrange} f(x)-P(x) = \frac{f^{[n+1]}(\xi_x)}{(n+1)!} \prod_{j=0}^n(x-x_j) \end{equation} \end{theorem} Ainsi l'erreur commise d\'epend d'une majoration de la taille de la d\'eriv\'ee $n+1$-i\`eme sur l'intervalle, mais aussi de la disposition des points $x_j$ par rapport \`a $x$. Par exemple si les points $x_j$ sont \'equidistribu\'es, le terme $|\prod_{j=0}^n(x-x_j)|$ sera plus grand pr\`es du bord de $I$ qu'au centre de $I$. Preuve du th\'eor\`eme~: Si $x$ est l'un des $x_j$ l'égalité est vraie. Soit \[ 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 \prod_{j=0}^n(t-x_j) \] elle s'annule en $x_j$ pour $j$ variant de 0 à $n$ ainsi qu'en $x$ suite au choix de la constante $C$, donc $g$ s'annule au moins $n+2$ fois sur l'intervalle contenant les $x_j$ et $x$, donc $g'$ s'annule au moins $n+1$ fois sur ce même intervalle, donc $g'{'}$ s'annule au moins $n$ fois, etc. et finalement $g^{[n+1]}$ s'annule une fois au moins sur cet intervalle. Or \[ g^{[n+1]} = f^{[n+1]} - C (n+1)!\] car $P$ est de degré inférieur ou égal à $n$ et $ \prod_{j=0}^n(x-x_j) - x^{n+1}$ est de degré inférieur ou égal à $n$. Donc il existe bien un réel $\xi_x$ dans l'intervalle contenant les $x_j$ et $x$ tel que \[ C=\frac{f^{[n+1]}(\xi_x)}{(n+1)!} \] {\bf Calcul efficace du polyn\^ome de Lagrange.}\\ Avec la m\'ethode de calcul pr\'ec\'edent, on remarque que le polyn\^ome de Lagrange peut s'\'ecrire \`a la Horner sous la forme~: \begin{eqnarray*} 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{eqnarray*} ce qui permet de le calculer rapidement une fois les $\alpha_i$ connus. On observe que \[ \alpha_0=f(x_0), \quad \alpha_1=\frac{f(x_1)-f(x_0)}{x_1-x_0} \] On va voir que les $\alpha_k$ peuvent aussi se mettre sous forme d'une diff\'erence. On d\'efinit les diff\'erences divis\'ees d'ordre $n$ par r\'ecurrence \[ 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 $\alpha_k=f[x_0,...,x_k]$. C'est vrai au rang 0, il suffit donc de le montrer au rang $k+1$ en l'admettant au rang $k$. Pour cela on observe qu'on peut construire le polyn\^ome d'interpolation en $x_0,...,x_{k+1}$ \`a partir des polyn\^omes d'interpolation $P_k$ en $x_0,...,x_k$ et $Q_k$ en $x_1,...,x_{k+1}$ par la formule~: \[ 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\'erifie que $P_{k+1}(x_i)=f(x_i)$ pour $i\in [1,k]$ car $P_k(x_i)=f(x_i)=Q_k(x_i)$, et pour $i=0$ et $i=k+1$, on a aussi $P_{k+1}(x_0)=f(x_0)$ et $P_{k+1}(x_{k+1})=f(x_{k+1})$. Or $\alpha_{k+1}$ est le coefficient dominant de $P_{k+1}$ donc c'est la diff\'erence du coefficient dominant de $Q_k$ et de $P_k$ divis\'ee par $x_{k+1}-x_0$, c'est-\`a-dire la d\'efinition de $f[x_0,...,x_{k+1}]$ en fonction de $f[x_1,...,x_{k+1}]$ et $f[x_0,...,x_{k}]$. Exemple~: on reprend $P(0)=1, P(1)=2, P(2)=1$. On a \[ \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)=\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_i$, en observant que le coefficient dominant de $P$ ne d\'epend pas de cet ordre, on en d\'eduit que $f[x_0,...,x_k]$ est ind\'ependant de l'ordre des $x_i$, on peut donc \`a partir du tableau ci-dessus \'ecrire $P$ par exemple avec l'ordre 2,1,0, sous la forme \[ P(x)=1+(x-2)(-1+(x-1)(-1))=1+(x-2)(-x)\] %TP: Approximation de Padé pour exp. %Programmation de la recherche de toutes les racines %d'un polynôme. Bassins d'attractions. \section{Intégration numérique} % équations différentielles(?) \index{integration} \index{quadrature} Les fractions rationnelles admettent une primitive que l'on calcule en décomposant la fraction avec Bézout comme expliqué précédemment. Mais elles font figure d'exceptions, la plupart des fonctions n'admettent pas de primitives qui s'expriment à l'aide des fonctions usuelles. Pour calculer une intégrale,on revient donc à la définition d'aire sous la courbe, aire que l'on approche, en utilisant par exemple un polynome de Lagrange. Le principe est donc le suivant~: on découpe l'intervalle d'intégration en subdivisions $[a,b]=[a,a+h] + [a+h,a+2h]+...[a+(n-1)h,a+nh=b$, où $h=(b-a)/n$ est le pas de la subdivision, et sur chaque subdivision, on approche l'aire sous la courbe. \subsection{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(\alpha)$, $f(\beta)$ (rectangle à droite et gauche) ou $f((\alpha+\beta)/2)$ (point milieu), pour les trapèzes on utilise le segment reliant $[\alpha,f(\alpha)]$ à $[\beta,f(\beta)]$. Exemple~: calcul de la valeur approchée de $\int_0^1 t^3 dt$ (on en connait la valeur exacte $1/4=0.25$) par ces méthodes en subdivisant $[0,1]$ en 10 subdivisions (pas $h=1/10$), donc $\alpha=j/10$ et $\beta=(j+1)/10$ pour $j$ variant de 0 à 9. Pour les rectangles à gauche, on obtient sur une subdivision $f(\alpha)=(j/10)^3 $ que l'on multiplie par la longueur de la subdivision soit $h=1/10$~: \[ \frac{1}{10} \sum_{j=0}^9 (\frac{j}{10})^3 = \frac{81}{400} = 0.2025\] Pour les rectangles à droite, on obtient \[ \frac{1}{10} \sum_{j=1}^{10} (\frac{j}{10})^3 = \frac{121}{400} = 0.3025\] Pour le point milieu $f((\alpha+\beta)/2)=f((j/10+(j+1)/10)/2)=f(j/10+1/20)$ \[ \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 $x$, les verticales $y=\alpha$, $y=\beta$ et les points sur ces verticales d'ordonnées respectives $f(\alpha)$ et $f(\beta)$ vaut \[ h \frac{f(\alpha)+f(\beta)}{2}\] donc \[ \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~: \index{rectangle} \index{trapeze} \index{point milieu} \begin{eqnarray} \mbox{rectangle gauche} & = & h \sum_{j=0}^{n-1} f(a+jh) \label{eq:rectangle_gauche}\\ \mbox{rectangle droit} & = & h \sum_{j=1}^{n} f(a+jh) \label{eq:rectangle_droit}\\ \mbox{point milieu} & = & h \sum_{j=0}^{n-1} f(a+jh+\frac{h}{2}) \label{eq:point_milieu}\\ \mbox{ trapezes } & = & h \left(\frac{f(a)+f(b)}{2}+\sum_{j=1}^{n-1} f(a+jh) \right) \label{eq:trapezes} \end{eqnarray} où $h=(b-a)/n$ est le pas de la subdivision, $n$ le nombre de subdivisions. On observe sur l'exemple que le point milieu et les trapèzes donnent une bien meilleure précision que les rectangles. Plus généralement, la précision de l'approximation n'est pas la même selon le choix de méthode. Ainsi pour les rectangles à gauche (le résultat est le même à droite), si $f$ est continument dérivable, de dérivée majorée par une constante $M_1$ sur $[a,b]$, en faisant un développement de Taylor de $f$ en $\alpha$, on obtient \[ |\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=3$, l'erreur est donc majorée par $0.015$ sur une subdivision, donc par $0.15$ sur les 10 subdivisions. Pour le point milieu, on fait le développement en $(\alpha+\beta)/2$ à l'ordre 2, en supposant que $f$ est deux fois continument dérivable~: \begin{eqnarray*} |\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{eqnarray*} Dans l'exemple, on a $M_2=6$, donc l'erreur sur une subdivision est majorée par $0.25e-3$, donc sur 10 subdivisions par $0.25e-2=0.0025$. Pour les trapèzes, la fonction $g$ dont le graphe est le segment reliant $[\alpha,f(\alpha)]$ à $[\beta,f(\beta)]$ est $f(\alpha)+(t-\alpha)/(\beta-\alpha)f(\beta)$, c'est en fait un polynome de Lagrange, si $f$ est deux fois continument dérivable, on peut donc majorer la différence entre $f$ et $g$ en utilisant (\ref{eq:lagrange}), on intègre la valeur absolue ce qui donne \[ |\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} \] où $M_2$ est un majorant de $|f'{'}|$ sur $[a,b]$. Lorsqu'on calcule l'intégrale sur $[a,b]$ par une de ces méthodes, on fait la somme sur $n=(b-a)/h$ subdivisions de longueur $\beta-\alpha=h$, on obtient donc une majoration de l'erreur commise sur l'intégrale~: \begin{itemize} \item pour les rectangles à droite ou gauche $nM_1h^2/2=M_1 h (b-a)/2$ \item pour le point milieu $M_2 h^2 (b-a)/24$ \item pour les trapèzes $M_2h^2 (b-a)/12$. \end{itemize} Lorsque $h$ tend vers 0, l'erreur tend vers 0, mais pas à la même vitesse, plus rapidement pour les trapèzes et le point milieu que pour les rectangles. Plus on approche précisément la fonction sur une subdivision, plus la puissance de $h$ va être grande, plus la convergence sera rapide lorsque $h$ sera petit, avec toutefois une contrainte fixée par la valeur de $M_k$, borne sur la dérivée $k$-ième de $f$ (plus $k$ est grand, plus $M_k$ est grand en général). Nous allons voir dans la suite comment se comporte cette puissance de $h$ en fonction de la facon dont on approche $f$. \subsection{Ordre d'une méthode} \index{ordre} On appelle méthode d'intégration l'écriture d'une approximation de l'intégrale sur une subdivision sous la forme \[ \int_{\alpha}^{\beta} f(t) dt \approx I(f)=\sum_{j=1}^k w_j f(y_j) \] où les $y_j$ sont dans l'intervalle $[\alpha,\beta]$, par exemple équirépartis sur $[\alpha,\beta]$. On utilise aussi la d\'efinition~: \[ \int_{\alpha}^{\beta} f(t) dt \approx I(f)= (\beta-\alpha)\sum_{j=1}^k \tilde{w}_j f(y_j) \] On prend toujours $\sum_j w_j=\beta-\alpha$ (ou $\sum_j \tilde{w}_j=1$) pour que la m\'ethode donne le r\'esultat exact si la fonction est constante. On dit qu'une méthode d'intégration est d'ordre $n$ si il y a égalité ci-dessus pour tous les polynômes de degré inférieur ou \'egal à $n$ et non \'egalit\'e pour un polyn\^ome de degr\'e $n+1$. Par exemple, les rectangles à droite et gauche sont d'ordre 0, le point milieu et les trap\`ezes sont d'ordre 1. Plus g\'en\'erallement, si on approche $f$ par son polyn\^ome d'interpolation de Lagrange en $n+1$ points (donc par un polyn\^ome de degr\'e inf\'erieur ou \'egal \`a $n$), on obtient une m\'ethode d'int\'egration d'ordre au moins $n$. Si une méthode est d'ordre $n$ avec des $w_j\geq 0$ et si $f$ est $n+1$ fois continument dérivable, alors sur une subdivision, on a~: \begin{equation} \label{eq:majoration_err} |\int_{\alpha}^{\beta} f-I(f)| \leq M_{n+1} \frac{(\beta-\alpha)^{n+2}}{(n+1)!} (\frac{1}{n+2}+1) \end{equation} En effet, on fait le développement de Taylor de $f$ par exemple en $\alpha$ à l'ordre $n$ \begin{eqnarray*} 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{eqnarray*} Donc \[ |\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, \begin{eqnarray*} |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{eqnarray*} Donc comme la méthode est exacte pour $T_n(f)$, on en déduit que \begin{eqnarray*} |\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{eqnarray*} Si les $w_j\geq 0$, alors $\sum_{j=1}^k |w_j|=\sum_{j=1}^k w_j=\beta-\alpha$ et on obtient finalement (\ref{eq:majoration_err}) On remarque qu'on peut am\'eliorer la valeur de la constante en faisant tous les d\'eveloppement de Taylor en $(\alpha+\beta)/2$ au lieu de $\alpha$, Après sommation sur les $n$ subdivisions, on obtient que~: \begin{theorem} Pour une méthode d'ordre $n$ à coefficients positifs et une fonction $f$ $n+1$ fois continument d\'erivable~ \[|\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)\] \end{theorem} On observe que cette majoration a la bonne puissance de $h$ sur les exemples déja traités, mais pas forcément le meilleur coefficient possible, parce que nous avons trait\'e le cas g\'en\'eral d'une m\'ethode d'ordre $n$. \subsection{Simpson} \index{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 $\alpha,(\alpha+\beta)/2,\beta$. On calcule l'int\'egrale par exemple avec un logiciel de calcul formel, avec Xcas~: \begin{center} \verb|factor(int(lagrange([a,(a+b)/2,b],[fa,fm,fb]),x=a..b))| \end{center} qui donne la formule sur une subdivision \[ I(f) = \frac{h}{6} (f(\alpha)+4f(\frac{\alpha+\beta}{2}) + f(\beta)) \] et sur $[a,b]$~: \begin{equation} \label{eq:simpson} 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) \end{equation} Si on intègre $t^3$ sur $[0,1]$ en 1 subdivision par cette méthode, on obtient \[ \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 $f$ polynome de degré inférieur ou égal à 2 puisque l'approximation de Lagrange de $f$ est alors égale à $f$. On en déduit que la méthode de Simpson est d'ordre 3 (pas plus car la méthode de Simpson appliquée à l'intégrale de $t^4$ sur $[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~: \[ |\int_a^b f - I(f)| \leq \frac{h^4}{2880} (b-a) M_4 \] Cette méthode nécessite $2n+1$ évaluations de $f$ (le calcul de $f$ est un point étant presque toujours l'opération la plus couteuse en temps d'une méthode de quadrature), au lieu de $n$ pour les rectangles et le point milieu et $n+1$ pour les trapèzes. Mais on a une majoration en $h^4$ au lieu de $h^2$ donc le ``rapport qualité-prix'' de la méthode de Simpson est meilleur, on l'utilise donc plutot que les m\'ethodes pr\'ec\'edentes sauf si $f$ n'a pas la r\'egularit\'e suffisante (ou si $M_4$ est trop grand). \subsection{Newton-Cotes} \index{Newton-Cotes} On peut généraliser l'idée précédente, découper la subdivision $[\alpha,\beta]$ en $n$ parts \'egales et utiliser le polynôme d'interpolation en ces $n+1$ points $x_0=\alpha, x_1, ..., x_n=\beta$. Ce sont les méthodes de Newton-Cotes, qui sont d'ordre $n$ au moins. Comme le polyn\^ome d'interpolation d\'epend lin\'eairement des ordonn\'ees, cette m\'ethode est bien de la forme~: \[ I(f)=(\beta-\alpha)\sum_{j=0}^n \tilde{w}_j f(x_j)\] De plus les $\tilde{w}_j$ sont universels (ils ne d\'ependent pas de la subdivision), parce qu'on peut faire le changement de variables $x=\alpha+t(\beta-\alpha)$ dans l'int\'egrale et le polyn\^ome d'interpolation et donc se ramener \`a $[0,1]$. Exemple~: on prend le polyn\^ome d'interpolation en 5 points \'equidistribu\'es sur une subdivision $[a,b]$ (m\'ethode de Boole). Pour calculer les $\tilde{w}_j$, on se ram\`ene \`a $[0,1]$, puis on tape \begin{center} \verb|int(lagrange(seq(j/4,j,0,4),[f0,f1,f2,f3,f4]),x=0..1)| \end{center} et on lit les coefficients de \verb|f0| \`a \verb|f4| qui sont les $\tilde{w}_0$ \`a $\tilde{w}_4$: 7/90, 32/90, 12/90, 32/90, 7/90. La m\'ethode est d'ordre au moins 4 par construction, mais on v\'erifie qu'elle est en fait d'ordre 5 (exercice), la majoration de l'erreur d'une m\'ethode d'ordre 5 est \[ |\int_a^b f -I(f)| \leq \frac{M_6}{2^6 6!}(1+\frac{1}{7}) h^6 (b-a) \] elle peut \^etre am\'elior\'ee pour cette m\'ethode pr\'ecise en \[ |\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 $n\geq 8$, les $w_j$ ne sont pas tous positifs, et d'autre part, parce que la constante $M_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]$, puis 2, puis $2^2$, ..., et on élimine les puissances de $h$ du reste $\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 $h$ ne contient que des puissances paires de $h$). De plus, on peut être amené à faire varier le pas $h$ en fonction de la plus ou moins grande régularité de la fonction. \subsection{En r\'esum\'e} Int\'egration sur $[a,b]$, $h$ pas d'une subdivision, $M_k$ majorant de la d\'eriv\'ee $k$-i\`eme de la fonction sur $[a,b]$\\ \begin{tabular}{|c|c|c|c|c|} \hline & formule & Lagrange degr\'e & ordre & erreur \\ \hline rectangles & (\ref{eq:rectangle_gauche}), (\ref{eq:rectangle_droit}) & 0 & 0 & $M_1 h (b-a)/2$ \\ point milieu & (\ref{eq:point_milieu}) & 0 & 1 & $M_2 h^2 (b-a)/24$ \\ trap\`ezes & (\ref{eq:trapezes}) & 1 & 1 & $M_2 h^2 (b-a)/12$ \\ Simpson & (\ref{eq:simpson})& 2 & 3 & $M_4 h^4 (b-a)/2880$ \\ \hline \end{tabular} % Romberg, Gauss. %Comparaison série/intégrale. %Séries entières. %Point milieu, RK. %Séries entières. \section{Algèbre linéaire} \subsection{Le pivot de Gauss}\index{pivot} \index{Gauss} Cet algorithme permet de créer des zéros en effectuant des manipulations réversibles sur les lignes d'une matrice. Ces lignes peuvent représenter les coefficients d'un système linéaire, on obtient alors un système linéaire équivalent, ou les coordonnées d'un système de vecteur, on obtient alors les coordonnées d'un système de vecteur engendrant le même sous-espace vectoriel. On peut également représenter 2 matrices $A$ et $B$ reliés par une relation $Ax=B$, cette relation reste alors vraie au cours et donc après la réduction. \subsubsection{L'algorithme} L'algorithme est le suivant: \begin{enumerate} \item on initialise $c=1$ et $l=1$, $c$ désigne le numéro de colonne $c$ à réduire, et $l$ le numéro de ligne à partir duquel on cherche un ``pivot'' (au début $l$ et $c$ valent donc 1, en général les 2 augmentent de 1 à chaque itération) \item Si $c$ ou $l$ est plus grand que le nombre de colonnes ou de lignes on arrête. \item Si la colonne $c$ n'a que des coefficients nuls à partir de la ligne $l$, on incrémente le numéro de colonne $c$ de 1 et on passe à l'étape 2. Sinon, on cherche la ligne dont le coefficient est en valeur absolue le plus grand possible (en calcul approché) ou le plus simple possible (en calcul exact), on échange cette ligne avec la ligne $l$. Puis on effectue pour toutes les lignes sauf $l$ ou pour toutes les lignes à partir de $l+1$ (selon qu'il s'agit d'une réduction de Gauss complète ou d'une réduction de Gauss sous-diagonale) la manipulation réversible \[ L_j \leftarrow L_j -\frac{a_{jc}}{a_{lc}} L_l \] On incrémente $c$ et $l$ de 1 et on passe à l'étape 2. \end{enumerate} \subsubsection{Efficacité de l'algorithme} Si la matrice possède $L$ lignes et $C$ colonnes, le nombre maximal d'opérations pour réduire une ligne est $C$ divisions, $C$ multiplications, $C$ soustractions, donc $3C$ opérations arithmétiques de base. Il y a $L-1$ lignes à réduire à chaque étape et min$(L,C)$ étapes à effectuer, on en déduit que le nombre maximal d'opérations pour réduire une matrice est $3LC$min$(L,C)$. Pour une matrice carrée de taille $n$, cela fait $3n^3$ opérations. \subsubsection{Erreurs d'arrondis du pivot de Gauss} \index{erreur} Comme $|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^n$ (où $n$ est le nombre d'étapes de réduction, inférieur au plus petit du nombre de lignes et de colonnes). Ceci signifie qu'avec la précision d'un double, on peut au pire perdre toute précision pour des matrices pas si grandes que ça ($n=52$). Heureusement, il semble qu'en pratique, l'erreur absolue ne soit que très rarement multipliée par un facteur supérieur à 10. Par 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. {\bf Exemple}\\ Soit \`a r\'esoudre le syst\`eme lin\'eaire \[ \epsilon x + 1.0 y = 1.0 , \quad x + 2.0 y = 3.0 \] avec $\epsilon =2^{-54}$ (pour une machine utilisant des doubles pour les calculs en flottant, plus g\'en\'erallement on choisira $\epsilon$ tel que $(1.0+3\epsilon)-1.0$ soit indistinguable de 0.0).\\ Si on r\'esoud le syst\`eme exactement, on obtient $x=1/(1-2\epsilon)$ (environ 1) et $y=(1-3\epsilon)/(1-2\epsilon)$ (environ 1). Supposons que l'on n'utilise pas la strat\'egie du pivot partiel, on prend alors comme pivot $\epsilon$, donc on effectue la manipulation de ligne $L_2 \leftarrow L_2 - 1/\epsilon L_1$ ce qui donne comme 2\`eme \'equation $(2.0-1.0/\epsilon)y=3.0-1.0/\epsilon$. Comme les calculs sont num\'eriques, et \`a cause des erreurs d'arrondis, cette 2\`eme \'equation sera remplac\'ee par $(-1.0/\epsilon)y=-1.0/\epsilon$ d'o\`u $y=1.0$, qui sera remplac\'e dans la 1\`ere \'equation, donnant $\epsilon x = 1.0-1.0y=0.0$ donc $x=0.0$.\\ Inversement, si on utilise la strat\'egie du pivot partiel, alors on doit \'echanger les 2 \'equations $L_2'=L_1$ et $L_1'=L_2$ puis on effectue $L_2 \leftarrow L_2' - \epsilon L_1'$, ce qui donne $(1.0-2.0\epsilon) y = 1.0 - 3.0 \epsilon $, remplac\'ee en raison des erreurs d'arrondi par $1.0*y=1.0$ donc $y=1.0$, puis on remplace $y$ dans $L_1'$ ce qui donne $x=3.0-2.0y=1.0$.\\ On observe dans les deux cas que la valeur de $y$ est proche de la valeur exacte, mais la valeur de $x$ dans le premier cas est grossi\`erement eloign\'ee 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 $A$ 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). \subsection{Applications de Gauss} \subsubsection{Base d'un sous-espace} On réduit la matrice des vecteurs écrits en ligne, puis on prend les lignes non nulles, dont les vecteurs forment une base du sous-espace vectoriel engendré par les vecteurs du départ. Exemple~: base du sous-espac engendré par $(1,2,3), (4,5,6), (7,8,9)$. On réduit la matrice, la 3ème ligne est nulle donc on ne garde que les 2 premières lignes $(1,0,-1), (0,1,2)$ (remarque: une réduction sous-diagonale aurait suffi). \subsubsection{Déterminant} \index{determinant} 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. \subsubsection{Réduction sous forme échelonnée (rref)} \index{rref} \index{reduction} 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=b$ où $x$ est la solution de l'équation, et à la fin de la réduction $A=I$. Par exemple pour résoudre le système \[ \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 \verb|[[1,2,3,5],[4,5,6,0],[7,8,0,1]]|, ce qui donne \verb|[[1,0,0,-9],[0,1,0,8],[0,0,1,-2/3]]|, d'où la solution $x=-9, y=8, z=-2/3$. \subsubsection{Inverse} \index{inverse} On accolle la matrice identité à droite de la matrice à inverser. On effectue la réduction sous forme échelonnée, on doit obtenir à droite l'identité si la matrice est inversible, on a alors à gauche la matrice inverse. La relation conservée est en effet $A x=B$ où $x$ est l'inverse de la matrice de départ, et en fin de réduction $A=I$. Par exemple, pour calculer l'inverse de \verb|[[1,2,3],[4,5,6],[7,8,0]]|, on réduit avec \verb|rref| \verb|[[1,2,3,1,0,0],[4,5,6,0,1,0],[7,8,0,0,0,1]]|. \subsubsection{Noyau} \index{noyau} \index{ker} 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 \verb|[[1,2,3,4],[1,2,7,12]]|, on réduit la matrice avec rref, ce qui donne \verb|[[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 \verb|[[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 \verb|[[1,2,0,-2],[0,0,0,0],[0,0,1,2],[0,0,0,0]]|. La 2ème colonne donne le premier vecteur de la base du noyau, $(2,-1,0,0)$, la 4ème colonne donne le deuxième vecteur $(-2,0,2,-1)$, on vérifie aisément que ces 2 vecteurs forment une famille libre du noyau, donc une base car la dimension du noyau est égale à 2 (dimension de l'espace de départ moins le rang de la matrice, c'est-à-dire le nombre de lignes non nulles de la matrice réduite). \subsubsection{La méthode de factorisation $LU$} \index{LU} Nous ne la développons pas à ce niveau, elle permet d'écrire une matrice $A$ 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. \subsection{Réduction exacte des endomorphismes} \index{diagonalisation} On calcule le polynome caractéristique ou le polynome minimal, on le factorise, et on calcule ensuite le noyau de $A-\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. \subsubsection{Polynome caractéristique} \index{polynome caracteristique} On peut le calculer en développant le déterminant det$A-\lambda I$, mais il est plus efficace de le calculer par interpolation. Soit $A$ une matrice carrée de taille $n$, on sait que son polynome caractéristique est un polynome de degré $n$, il suffit de connaitre sa valeur en $n+1$ points distincts, on calcule donc $n+1$ déterminants det$A-\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 \verb|[[1,-1],[2,4]]|, calculer det$(A-\lambda I$) en $\lambda=0,1,2$ puis le polynome d'interpolation, vérifier que c'est bien le polynome caractéristique. Il faut effectuer $n+1$ calculs de déterminants, ce qui nécessite $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)$ opérations). \subsubsection{Polynome minimal} \index{polynome minimal} \begin{defn} Le polynome minimal d'une matrice $A$ est un polynôme $M$ de degré minimal tel que $M(A)=0$ et de coefficient dominant égal à 1. Un tel polynome divise tous les polynomes tels que $P(A)=0$, il divise le polynome caractéristique de $A$ et il a les mêmes racines que le polynome caractéristique. \end{defn} Preuve:\\ D'abord $M$ divise tous les polynomes tels que $P(A)=0$, car si $R$ désigne le reste de la division de $P$ par $M$ alors $R(A)=(P-QM)(A)=P(A)-Q(A)M(A)=0$, donc $R$ est nul car son degré est plus petit que celui de $M$. En particulier le polynome minimal divise le polynome caractéristique $C$, car $C(A)=0$ (on peut montrer que $C(A)=0$ en faisant le produit de la matrice $A-\lambda I$ par sa comatrice, on obtient le déterminant fois l'identité, soit $C(\lambda)I$. Comme $C(\lambda)I-C(A)$ peut se factoriser par $\lambda I -A$ en appliquant (\ref{eq:factor}) à chaque monome de $C$, on en déduit que $C(A)$ se factorise par $\lambda I -A$, donc $C(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-\lambda I$ n'est pas inversible. Or $M(A)-M(\lambda)I$ se factorise par $A-\lambda I$ car \begin{equation} \label{eq:factor} A^k-\lambda^k I=(A-\lambda I) \sum_{j=0}^{k-1} \lambda^{k-1-j} A^j \end{equation} donc $M(A)-M(\lambda)I$ ne peut pas être inversible. Comme $M(A)=0$ on en déduit que $M(\lambda)I$ n'est pas inversible donc $M(\lambda)=0$, $\lambda$ est une racine de $M$. Donc si le polynome caractéristique n'a pas de racines multiples, il est égal au polynome minimal. Pour calculer $M$, on peut chercher une relation de degré minimal entre les puissances de $A$, en les voyant comme des vecteurs à $n^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 $A$ (de 0 à $n$), 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^2$ lignes et $n+1$ colonnes. Il existe une autre méthode moins couteuse et qui marche presque toujours. Elle consiste à calculer le polynome minimal de $A$ par rapport à un vecteur $v$ c'est-à-dire le polynome de degré minimal (et coefficient dominant 1) tel que $M_v(A)v=0$. Comme $M(A)=0$, on a $M(A)v=0$, donc $M_v$ divise $M$ qui divise le polynome caractéristique. Si par chance, on trouve que $M_v$ est de degré $n$, alors $M_v$ sera égal à $M$ et au polynome caractéristique. On fait donc le calcul du noyau de l'application linéaire dont les colonnes sont les $A^jv$ pour $j$ variant de 0 à $n$. Si l'on trouve un espace de dimension 1, alors $M_v$ est de degré $n$ 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 $n$, on peut essayer un ou quelques autres vecteurs, et faire le PPCM des polynomes minimaux obtenus. Si on obtient un polynome de degré $n$ on conclut, sinon on peut tester si ce polynome évalué en $A$ est nul, ce sera alors le polynome minimal. Exemple, on reprend la matrice \verb|[[1,-1],[2,4]]|, et comme vecteur aléatoire $v=(1,0)$, on a $Av=(1,-1)$ et $A(Av)=(-1,-5)$. On calcule donc le noyau de la matrice \verb|[[1,1,-1],[0,-1,-5]]| (on écrit en colonnes $v, Av, A^2v$), on trouve que $(-6,5,-1)$ engendre le noyau, donc le polynome minimal relatif au vecteur $v$ est (au signe près) $-6+5x-x^2$. Comme il est de degré maximal 2, c'est le polynome minimal et caractéristique. \subsection{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-\lambda I$), lorsqu'on veut calculer quelques valeurs propres on préfère utiliser des méthodes itératives directement sur $A$ ce qui évite la propagation des erreurs d'arrondi. \subsubsection{Méthode de la puissance} \index{puissance} Elle permet de déterminer la plus grande valeur propre en valeur absolue d'une matrice diagonalisable lorsque celle-ci est unique. Supposons en effet que les valeurs propres de $A$ soient $x_1,...,x_n$ avec $|x_1| \leq |x_2| \leq ... \leq |x_{n-1}| < |x_n|$ et soient $e_1,...,e_n$ une base de vecteurs propres correspondants. On choisit un vecteur aléatoire $v$ et on calcule la suite $v_n=Av_{n-1}=A^n v $. Si $v$ a pour coordonnées $V_1,...,V_n)$ dans la base propre, alors \[ 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_n$ est l'unique valeur propre de module maximal entraine alors que $\lim_{n \rightarrow +\infty} w_n = V_n e_n$ puisque la suite géométrique de raison $x_j/x_n$ converge vers 0. Autrement dit, si $V_n\neq 0$ (ce qui a une probabilité 1 d'être vrai pour un vecteur aléatoire), $v_n$ est équivalent à $V_n x_n^n e_n$. Lorsque $n$ est grand, $v_n$ est presque colinéaire au vecteur propre $e_n$ (que l'on peut prendre comme $v_n$ divisé par sa norme), ce que l'on détecte en testant si $v_{n+1}$ et $v_n$ sont presques colinéaires, et de plus le facteur de colinéarité entre $v_{n+1}$ et $v_n$ est presque $x_n$, la valeur propre de module maximal. Exercice~: tester la convergence de $A^nv$ vers l'espace propre associé à $\lambda=3$ pour la matrice \verb|[[1,-1],[2,4]]| et le vecteur $(1,0)$. Lorsqu'on applique cette méthode a une matrice réelle, il peut arriver quíl y ait deux valeurs propres conjuguées de module maximal. Le même type de raisonnement montre que pour $n$ grand, $v_{n+2}$ est presque colinéaire à l'espace engendré par $v_n$ et $v_{n+1}$, la relation $ 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. \subsubsection{Itérations inverses} \index{iterations inverses} La méthode précédente permet de calculer la valeur propre de module maximal d'une matrice. Pour trouver une valeur propre proche d'une quantité donnée $x$, on peut appliquer la méthode précédente à la matrice $(A-xI)^{-1}$. En effet, les valeurs propres de cette matrice sont les $(x_i-x)^{-1}$ dont la norme est maximale lorsqu'on se rapproche de $x_i$. \subsubsection{Elimination des valeurs propres trouvées} Si la matrice $A$ est symétrique, et si $e_n$ est un vecteur propre normé écrit en colonne, on peut considérer la matrice $B=A-x_n e_n e_n^t$ qui possède les mêmes valeurs propres et mêmes vecteurs propres que $A$ avec même multiplicité, sauf $x_n$ qui est remplacé par 0. En effet les espaces propres de $A$ sont orthogonaux entre eux, donc \[ 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 $A$ n'est pas symétrique, on peut utiliser une technique analogue lorsque 0 n'est pas valeur propre de $A$ (on peut s'y ramener en ajoutant à $A$ un multiple de l'identité). En effet on peut construire un vecteur propre de $B$ pour une valeur propre $x_k \neq 0$ à partir d'un vecteur propre de $B$, en cherchant $y$ tel que tel que \[ B(e_k-ye_n)=x_k(e_k-ye_n) \] On obtient pour le membre de gauche~: \[ 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_k e_k - y x_k e_n \] d'où l'équation \[ 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_k$ est proche de 0, car les vecteurs propres se rapprochent alors tous de $e_n$. %TP: cf tpmias2 %TP: a faire \section{Quelques références} \begin{itemize} \item Analyse numérique et équations différentielles, Demailly J.-P., Presses Universitaires de Grenoble, 1996 \item The Art of Computer Programming, Vol. 2: Seminumerical algorithms, Knuth D., Addison-Wesley, 1998 \item Mathématiques concrètes, illustrées par la TI-92 et la TI-89. Lemberg H, et Ferrard J.-M., Springer, 1998 \item Maths et Maple, J.M. Ferrard, Dunod, 1998 \item Handbook of Mathematical Functions, Abramowitz and Stegun, disponible en ligne sur \\ \verb|http://www.math.sfu.ca/~cbm/aands/toc.htm| \item Arithm\'etique flottante, Rapport de l'INRIA de V. Lef\`evre et P. Zimmermann, t\'el\'echargeable sur\\ \verb|http://www.inria.fr/rrrt/rr-5105.html| \item Modern Computer Arithmetic, R. P. Brent, P. Zimmermann, t\'el\'echargeable sur\\ \verb|http://www.loria.fr/~zimmerma/mca/pub226.html| \item Matrix computations, Golub and Loa, Hopkins University Press, 1989 \item Gantmacher \end{itemize} \pagebreak \printindex \pagebreak \appendix \section{La moyenne arithm\'etico-g\'eom\'etrique.} \subsection{D\'efinition et convergence} Soient $a$ et $b$ deux r\'eels positifs, on d\'efinit les 2 suites \begin{equation} \label{eq:agm} u_0=a, v_0=b, \quad u_{n+1}=\frac{u_n+v_n}{2}, v_{n+1}=\sqrt{u_nv_n} \end{equation} On va montrer que ces 2 suites sont adjacentes et convergent donc vers une limite commune not\'ee $M(a,b)$ et il se trouve que la convergence est tr\`es rapide, en raison de l'identit\'e~: \begin{equation} \label{eq:agm1} 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 \end{equation} la convergence est quadratique. On suppose dans la suite que $a\geq b$ sans changer la généralité puisque échanger $a$ et $b$ ne change pas la valeur de $u_n$ et $v_n$ pour $n>0$. On a alors $u_n \geq v_n$ (d'après (\ref{eq:agm1}) pour $n>0$) et $u_{n+1} \leq u_n$ car \[ u_{n+1}-u_n=\frac{1}{2}(v_n-u_{n}) \leq 0\] et $v_{n+1}=\sqrt{u_nv_n} \geq \sqrt{v_nv_n}=v_n$. Donc $(u_n)$ est décroissante minorée (par $v_0$), $(v_n)$ est croissante majorée (par $u_0$), ces 2 suites sont convergentes et comme $u_{n+1}=\frac{u_n+v_n}{2}$, elles convergent vers la même limite $l$ qui d\'epend de $a$ et $b$ et que l'on note $M(a,b)$. On remarque aussi que $M(a,b)=bM(a/b,1)=aM(1,b/a)$. Précisons maintenant la vitesse de convergence lorsque $a \geq b>0$. On va commencer par estimer le nombre d'itérations nécessaires pour que $u_n$ et $v_n$ soient du même ordre de grandeur. Pour cela, on utilise la majoration \[ \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 \[ \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 $n \geq \frac{\ln( \ln(a/b)/m)}{\ln(2)}$ alors $\ln \frac{u_n}{v_n} \leq m$ (par exemple, on peut prendre $m=0.1$ pour avoir $u_n/v_n \in [1,e^{0.1}])$. Le nombre minimum d'itérations $n_0$ est proportionnel au log du log du rapport $a/b$. Ensuite on est ramené à étudier la convergence de la suite arithmético-géométrique de premiers termes $a=u_{n_0}$ et $b=v_{n_0}$ et même en tenant compte de $M(a,b)=aM(1,b/a)$ à $a=1$ et $b=v_n/u_n$ donc $0\leq a-b \leq 1-e^{-0.1}$. Alors l'équation (\ref{eq:agm1}) entraine \[ u_{n+1}-v_{n+1} \leq \frac{1}{8}(u_n-v_n)^2 \] puis (par récurrence) \[ 0 \leq u_n-v_n \leq \frac{1}{8^{2^n-1}}(a-b)^{2^n} \] Donc comme $M(a,b)$ est compris entre $v_n$ et $u_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/\epsilon))$. Typiquement dans la suite, on souhaitera calculer $M(1,b)$ avec $b$ de l'ordre de $2^{-n}$ en déterminant $n$ chiffres significatifs, il faudra alors $O(\ln(n))$ itérations pour se ramener à $M(1,b)$ avec $b\in [e^{-0.1},1]$ puis $O(\ln(n))$ itérations pour avoir la limite avec $n$ chiffres significatifs. {\bf Le cas complexe}\\ On suppose maintenant que $a, b \in \C$ avec $\Re(a)>0, \Re(b)>0$. On va voir que la suite arithmético-géométrique converge encore. \\ {\bf \'Etude de l'argument}\\ On voit aisément (par récurrence) que $\Re(u_n)>0$ ; de plus $\Re(v_n) > 0$ car par définition de la racine carrée $\Re(v_n)\geq 0$ et est de plus non nul car le produit de deux complexes d'arguments dans $]-\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)$ se trouve dans l'intervalle de bornes $\arg(u_n)$ et $\arg(v_n)$ et que $\arg(v_{n+1})=\frac{1}{2}(\arg(u_n)+\arg(v_n))$ donc \[ | \arg(u_{n+1}-\arg(v_{n+1}) | \leq \frac{1}{2}|\arg(u_n)-\arg(v_n)| \] Après $n$ itérations, on a \[ |\arg(u_n)-\arg(v_n)| \leq \frac{\pi}{2^n} \] Après quelques itérations, $u_n$ et $v_n$ seront donc presque alignés. Faisons 4 itérations. On peut factoriser par exemple $v_n$ et on est ramené à l'étude de la suite de termes initiaux $a=u_n/v_n$ d'argument $\arg(u_n)-\arg(v_n)$ petit (inférieur en valeur absolue à $\pi/16$) et $b=1$. On suppose donc dans la suite que \[ |\arg(\frac{u_n}{v_n})| \leq \frac{\pi/16}{2^n} \] {\bf \'Etude du module}\\ On a~: \[ \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 $ \frac{u_{n}}{v_{n}}=\rho_n e^{i\theta_n}$, on a~: \begin{eqnarray*} |\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{eqnarray*} Si $\rho$ désigne le max de $\rho_n$ et $1/\rho_n$, on a alors la majoration \[ |\frac{u_{n+1}}{v_{n+1}}| \leq \frac{1}{2} \sqrt{ \rho + \rho + 2 \rho } = \sqrt{\rho} \] donc en prenant les logarithmes \begin{equation} \label{eq:lnrhon} \ln \rho_{n+1} \leq \frac{1}{2} \ln \rho=\frac{1}{2} |\ln \rho_n| \end{equation} On rappelle qu'on a la majoration \[ |\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 $\rho_{n+1}$ \begin{eqnarray*} \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{eqnarray*} en prenant les log et en minorant $\ln(1-x)$ par $-2x$ \[ \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 (\ref{eq:lnrhon}) \[ |\ln \rho_{n+1}| \leq \frac{1}{2} (|\ln \rho_n|+\frac{1}{2^{2n+3}} ) \] On en déduit \[ |\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)$ vers 0 est donc géométrique, donc $u_n$ et $v_n$ convergent quadratiquement. \subsection{Lien avec les int\'egrales elliptiques} Le calcul de la limite commune des suites $u_n$ et $v_n$ en fonction de $a$ et $b$ n'est pas trivial au premier abord. Il est reli\'e aux int\'egrales elliptiques, plus pr\'ecis\'ement on peut construire une int\'egrale d\'ependant de deux param\`etres $a$ et $b$ et qui est invariante par la transformation $u_n,v_n \rightarrow u_{n+1},v_{n+1}$ (\ref{eq:agm}) \[ I(a,b)=\int_{-\infty}^{+\infty} \frac{dt} {\sqrt{(a^2+t^2)(b^2+t^2)}} \] On a en effet \[ 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=\frac{1}{2} (t-\frac{ab}{t}), \quad t>0 \] o\`u $t \rightarrow u$ est une bijection croissante de $t\in]0,+\infty[$ vers $u \in ]-\infty,+\infty[$, donc \begin{eqnarray*} 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{eqnarray*} On note au passage que $I$ est définie si $a,b \in \C$ vérifient $\Re(a)>0, \Re(b)>0$, on peut montrer que la relation ci-dessus s'étend (par holomorphie). Lorsque $a=b=l$ (par exemple lorsqu'on est \`a la limite), le calcul de $I(l,l)$ est explicite \[ 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))=\frac{\pi}{M(a,b)}\] On peut transformer $I(a,b)$ en posant $t=bu$ \[ 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)$ ($du=(1+u^2) dx$) \[ 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)=\frac{\sin(x)^2}{1-\sin(x)^2}$ \[ 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\'efinit pour $m<1$ \[ K(m)=\int_0^{\frac{\pi}{2}} \frac{dx}{\sqrt{1-m \sin(x)^2}} \] alors on peut calculer $K$ en fonction de $I$, en posant $m=1-b^2/a^2$ soit $b^2/a^2=1-m$ \[ 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\`u l'on d\'eduit la valeur de l'int\'egrale elliptique en fonction de la moyenne arithm\'etico-g\'eom\'etrique~: \begin{equation} \label{eq:K} K(m)=\int_0^{\frac{\pi}{2}} \frac{dx}{\sqrt{1-m \sin(x)^2}}= \frac{\pi}{2M(1,\sqrt{1-m})} \end{equation} Dans l'autre sens, pour $x$ et $y$ positifs \[ 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)=\frac{\pi}{4} \frac{x+y}{ K\left( (\frac{x-y}{x+y}\right)^2 )}\] \subsection{Application~: calcul efficace du logarithme.} On peut utiliser la moyenne arithm\'etico-g\'eom\'etrique pour calculer le logarithme efficacement, pour cela on cherche le d\'eveloppement asymptotique de $K(m)$ lorsque $m$ tend vers 1. Plus pr\'ecis\'ement, on va poser $1-m=k^2$ avec $k \in ]0,1]$, donc \[ 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=\pi/2-x$, et \[ K(m)=\int_0^{\frac{\pi}{2}} \frac{dy}{\sqrt{\sin(y)^2+k^2 \cos(y)^2}}\] la singularit\'e de l'int\'egrale pour $k$ proche de 0 apparait lorsque $y$ est proche de 0. Si on effectue un d\'eveloppement de Taylor en $y=0$, on trouve \[ \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)$ \`a l'int\'egrale \[ 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=\frac{k}{\sqrt{1-k^2}} \sinh(t)\] ou directement avec Xcas, \begin{center} \verb|supposons(k>0 && k<1);|\\ \verb|J:=int(1/sqrt(k^2+(1-k^2)*y^2),y,0,pi/2)| \end{center} qui donne apr\`es r\'e\'ecriture~: \begin{equation} \label{eq:J} 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) \end{equation} et on peut calculer le d\'eveloppement asymptotique de $J$ en 0 \begin{center} \verb|series(J,k=0,5,1)| \end{center} qui renvoie~: \[ J =\ln\left(\frac{\pi}{k}\right) +O( \left(\frac{-1}{\ln(k)}\right)^5)\] on peut alors pr\'eciser ce d\'eveloppement par \begin{center} \verb|series(J+ln(k)-ln(pi),k=0,5,1)| \end{center} qui renvoie (apr\`es simplifications et où la notation $\tilde{O}$ peut contenir des logarithmes) \[ \left(\frac{1}{\pi^2} + \frac{\ln(\pi)-\ln(k)-1}{2}\right) k^{2} + \tilde{O}(k^4)\] donc \begin{equation} J=-\ln(k)+\ln(\pi)+\left(\frac{1}{\pi^2} + \frac{\ln(\pi)-\ln(k)-1}{2}\right) k^{2} + \tilde{O}(k^4) \end{equation} Examinons maintenant $K-J$, il n'a plus de singularit\'e en $y=0$, et il admet une limite lorsque $k \rightarrow 0$, obtenue en remplacant $k$ par 0 \[ (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\`u pour $K$ \[ K_{k \rightarrow 0} = \ln\left(\frac{4}{k}\right) + O( \left(\frac{-1}{\ln(k)}\right)^5)\] Pour pr\'eciser la partie du d\'eveloppement de $K$ en puissances de $k$, nous allons majorer $K-J-\ln(4/\pi)$, puis $J-\ln(\pi/k)$. Posons \[ A=\sin(y)^2+k^2 \cos(y)^2, \quad B=y^2+(1-y^2)k^2\] {\bf Majoration de $K-J-\ln(4/\pi)$}\\ L'int\'egrand de la diff\'erence $K-J-\ln(\frac{4}{\pi})$ est \begin{eqnarray} \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)} \label{eq:maj0} \\ &=& \frac{(y^2-\sin(y)^2)(1-k^2)}{\sqrt{A} \sqrt{B} (\sqrt{A}+\sqrt{B})} - \frac{y-\sin(y)}{y\sin(y)} \label{eq:maj1} \end{eqnarray} Soit \begin{equation} \label{eq:maj2} 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)} \end{equation} On décompose l'intégrale en 2 parties $[0,k]$ et $[k,\pi/2]$. Sur $[0,k]$ on utilise (\ref{eq:maj0}), on majore chaque terme séparément et on minore $A$ et $B$ par \[ 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 \begin{eqnarray*} | \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{eqnarray*} Sur $[k,\pi/2]$, on utilise (\ref{eq:maj2}) et on minore $A$ et $B$ par \[ 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 \[ | \int_k^{\frac{\pi}{2}} | \leq \int_k^{\frac{\pi}{2}} \frac{(y-\sin(y))|C|} {y \sin(y) (y+\sin(y))} , \] où~: \begin{eqnarray*} 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{eqnarray*} Donc \begin{eqnarray*} |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{eqnarray*} et \[ | \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 $y-\sin(y) \leq y^3/6$, donc \[ | \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 $A$ et $B$ par 1, \[ | \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 \begin{eqnarray*} \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{eqnarray*} Le deuxième morceau se majore en minorant $\sin(y)\geq (2y)/\pi$ \[ \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 \[ |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) \] où $J$ est donné en (\ref{eq:J}). {\bf Majoration de $J-ln(\pi/k)$}\\ On a \[ |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 $k\leq 1/2$, on a \[ \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 $\frac{1}{\sqrt{1-k^2}}$ par $\frac{2}{\sqrt{3}}$, l'argument du logarithme est inférieur à 1 et supérieur à \[ \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}) > 1-k^2 \] donc le logarithme en valeur absolue est inférieur à \[ 2 k^2 \] donc, pour $k\leq 1/2$, \[ |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/2$ \begin{equation} \label{eq:ln_agm0} |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) \end{equation} que l'on peut réécrire \begin{equation} \label{eq:ln_agm} |\frac{\pi}{2M(1,k)}-\ln\left(\frac{4}{k}\right) | \leq k^2(3.8-0.8\ln(k)) \end{equation} La formule (\ref{eq:ln_agm}) permet de calculer le logarithme d'un r\'eel positif avec (presque) $n$ bits lorsque $k \leq 2^{-n/2}$ (ce \`a quoi on peut toujours se ramener en calculant le logarithme d'une puissance $2^m$-i\`eme de $x$ ou le logarithme de $2^{m}x$, en calculant au pr\'ealable $\ln(2)$). Par exemple, prenons $k=2^{-27}$, on trouve (en 8 itérations) $M(1,2^-{27})=M_1=0.0781441403763$. On a, avec une erreur inférieure à $19 \times 2^{-54}=1.1\times 10^{-15}$ \[ M(1,2^-{27})=M_1=\frac{\pi}{2\ln(2^{29})}=\frac{\pi}{58\ln(2)}, \] On peut donc d\'eduire une valeur approch\'ee de $\pi $ si on connait la valeur approchée de $\ln(2)$ et réciproquement. Si on veut calculer les deux simultan\'ement, comme les relations entre $\ln$ et $\pi$ seront des \'equations homog\`enes, on est oblig\'e d'introduire une autre relation. Par exemple pour calculer une valeur approch\'ee de $\pi$ on calcule la diff\'erence $\ln(2^{29}+1)-\ln(2^{29})$ dont on connait le d\'eveloppement au premier ordre, et on applique la formule de la moyenne arithm\'etico-g\'eom\'etrique. Il faut faire attention \`a la perte de pr\'ecision lorsqu'on fait la diff\'erence des deux logarithmes qui sont tr\`es proches, ainsi on va perdre une trentaine de bits, il faut grosso modo calculer les moyennes arithm\'etico-g\'eom\'etrique avec 2 fois plus de chiffres significatifs. L'int\'er\^et de cet algorithme apparait lorsqu'on veut calculer le logarithme avec beaucoup de pr\'ecision, en raison de la convergence quadratique de la moyenne arithm\'etico-g\'eom\'etrique (qui est nettement meilleure que la convergence lin\'eaire pour les d\'eveloppements en s\'erie, 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, \`a partir du logarithme, ou les fonctions trigonom\'etriques inverses (en utilisant des complexes) et directes. On trouvera dans Brent-Zimmermann quelques consid\'erations permettant d'am\'eliorer les constantes dans les temps de calcul par rapport \`a cette m\'ethode (cela n\'ecessite d'introduire des fonctions sp\'eciales $\theta$) et d'autres formules pour calculer $\pi$. \end{document} %%% Local Variables: %%% mode: latex %%% TeX-master: t %%% End: