% Merci a Yannick Chevallier pour la version mathjax-enabled
\documentclass{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath,amsfonts,amssymb}
\usepackage[francais]{babel}
\usepackage{times}
\usepackage{xspace}
\usepackage{makeidx}
\usepackage{graphicx}
\usepackage{ifpdf}
\usepackage{a4wide}
\ifpdf
 \usepackage[colorlinks,pdftex]{hyperref}
\else
 \usepackage[ps2pdf,breaklinks=true,colorlinks=true,linkcolor=red,citecolor=green]{hyperref}
\fi
\newenvironment{preuve}{{\bf Preuve~:}}{}

%\newcommand{\bs}{\symbol{92}}
\newtheorem{prop}{Proposition}[section]
\newtheorem{thm}[prop]{Th\'{e}or\`{e}me}
\newtheorem{theorem}[prop]{Th\'{e}or\`{e}me}
\newtheorem{lem}[prop]{Lemme}
\newtheorem{coro}[prop]{Corollaire}
\newtheorem{defi}[prop]{D\'efinition}


%\theoremstyle{definition}

\newtheorem{defn}[prop]{D\'{e}finition}
\newtheorem{ex}[prop]{Exemple}
\newtheorem{exs}[prop]{Exemples}
\newtheorem{exe}[prop]{Exercice}
\newtheorem{rem}[prop]{Remarque}
\newtheorem{rems}[prop]{Remarques}
\newcommand{\trick}{\begin{flushleft}
\medskip
\end{flushleft}}
\makeindex

% \renewcommand{\cuttingunit}{subsection} (add HEVEA after % for hacha)
\input{giacfr.tex}
\ifhevea
\(
\newcommand{\R}{{\mathbb{R}}}
\newcommand{\mathbb{R}}{{\mathbb{R}}}
\newcommand{\C}{{\mathbb{C}}}
\newcommand{\mathbb{C}}{{\mathbb{C}}}
\newcommand{\Z}{{\mathbb{Z}}}
\newcommand{\N}{{\mathbb{N}}}
\newcommand{\Q}{{\mathbb{Q}}}
\newcommand{\tr}{\mbox{tr\,}}
\newcommand{\M}{\mathrm{M}}
\newcommand{\Mat}{\mathrm{Mat}}
\newcommand{\ds}{\displaystyle}
\)
\else
\newcommand{\R}{{\mathbb{R}}}
\newcommand{\rr}{{\mathbb{R}}}
\newcommand{\C}{{\mathbb{C}}}
\newcommand{\cc}{{\mathbb{C}}}
\newcommand{\Z}{{\mathbb{Z}}}
\newcommand{\N}{{\mathbb{N}}}
\newcommand{\Q}{{\mathbb{Q}}}
\newcommand{\tr}{\mbox{tr\,}}
\newcommand{\M}{\mathrm{M}}
\newcommand{\Mat}{\mathrm{Mat}}
\newcommand{\ds}{\displaystyle}
\fi
%HEVEA\htmlfoot{Retour \`a la page principale de \ahref{http://www-fourier.univ-grenoble-alpes.fr/\home{parisse}/giac_fr.html}{Giac/Xcas}.}
\giacmathjax
%HEVEA\renewcommand{\footertext}{}

\begin{document}
\title{Mat 406}
\author{Bernard.Parisse@univ-grenoble-alpes.fr}
\date{2022}

\begin{giacjshere}

\maketitle

\begin{center}
\includegraphics{mat406_qrcode.png}
\end{center}

N.B.~: la version HTML de ce document est interactive, vous pouvez
ex\'ecuter certains champs d'entr\'ee donn\'es en exemple
en les modifiant \'eventuellement, il suffit de
cliquer sur le bouton exe. Il est recommand\'e
d'utiliser Firefox ou un navigateur compatible supportant MathML.

\tableofcontents

\section{Pr\'esentation 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\'eances, 18h de TD en 12
s\'eances et 19.5h de TP en 13 s\'eances). 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 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 Interpolation polynomiale
\item Int\'egration approch\'ee.
\item Alg\`ebre lin\'eaire.
\end{enumerate}
L'ordre des deux premiers th\`emes sera invers\'e en cours.

L'\'evaluation se fait sur~:
\begin{itemize}
\item 1/4~: un DS \`a mi-semestre
\item 1/4~: certains compte-rendus de TP (\`a r\'ediger seul ou en binome), 
\item 1/2~: l'examen final
\end{itemize}
Les calculatrices 
%et les netbooks de taille d'\'ecran plus petits que 13 pouces 
sont autoris\'ees au DS et \`a l'examen final (pr\^et possible
pour le semestre).

\section{Repr\'esentation des nombres et autres donn\'ees, calcul exact/approch\'e}
R\'esum\'e:\\
Types de base~: entier machine, entier long, flottant machine et 
multipr\'ecision (Base 2, base 10, BCD). \\
Types compos\'es~: complexes, polynomes (repr\'esentation dense/creuse), 
symboles, listes (vecteurs, matrices), expressions, fonctions.\\
Erreur relative, erreur absolue, erreur d'arrondi, +/-, */\%\\
Algorithme, complexit\'e, exemple puissance modulaire, algorithme de Horner.

Les principaux ensembles de nombres en math\'ematiques sont les entiers
positifs $\mathbb{N}$ et relatifs $\mathbb{Z}$, les rationnels $\mathbb{Q}$, les r\'eels $\mathbb{R}$
et les complexes $\mathbb{C}$. Sur ordinateur, on peut repr\'esenter ces
nombres de mani\`ere exacte dans certains cas, approch\'ee dans d'autres.

\subsection{Représentation des entiers}

\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 \geq 0$.\\
Exemple~: \giacinput{iquorem(23,7)}

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 {\bf base} $b=10$\index{base}
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|0x19| car 25 divis\'e
par 16 donne quotient 1, reste 9\\
\giacinput{convert(25,base,16)}\\
En base 2, on trouverait \verb|0b11001| car $25=2^4+2^3+1$.\\
\giacinput{convert(25,base,2)}\\
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]$. 

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^{30}+2^{30}+2^{30}+2^{30}$
renverra 0. Ils sont utilisables avec tous les langages de
programmation traditionnels.

Les logiciels de calcul formel et certains logiciels
de programmation permettent de travailler avec
des entiers de taille beaucoup plus grande, ainsi qu'avec des
rationnels, permettant de faire du calcul exact, mais on paie
cette exactitude par un temps de calcul plus long, de plus
pas mal de m\'ethodes num\'eriques ne gagnent rien \`a
faire des calculs interm\'ediaires exacts. N\'eanmoins, l'utilisation
d'un logiciel de calcul formel permettra dans certains cas d'illustrer
certains ph\'enom\`enes dus au calcul approch\'e.

\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 {\bf mantisse}\index{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}\index{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'{\bf arrondi}\index{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, 
(puisque le test $x==y$ renvoie faux), on introduit les flottants 
{\bf d\'enormalis\'es} \index{d\'enormalis\'e}, 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\'eralement 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 {\bf BCD}\index{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 selon la norme IEEE-754, 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.5 - 5*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|

Voici un programme C++ affichant la 
\ahref{http://www-fourier.univ-grenoble-alpes.fr/\home{parisse}/mat406/format.cc}{repr\'esentation interne des flottants}

\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 arrondit 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} \index{erreur}
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.
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
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 et une plage d'exposants plus grande. Les calculs sont plus longs mais les erreurs
plus faibles. Attention, il s'agit toujours de calcul approch\'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. Dans ce cas, il est pertinent lorsqu'on
\'evalue $f(x)$ avec $x$ mal connu de calculer aussi $f'(x)$, en effet~:
$$ f(x(1+h))= f(x)+xh f'(x) + O(h^2)$$
l'erreur relative sur $f(x)$ est donc au premier ordre multipli\'ee par
$$ |\frac{xf'(x)}{f(x)}|$$
Par exemple, l'erreur relative sur $e^x$ est au premier ordre l'erreur relative sur $x$
multipli\'ee par $|x|$.
\giacinputbigmath{a:=convert(100,interval); (right(a)-left(a))/evalf(a,12)}\\
\giacinputbigmath{b:=exp(a); (right(b)-left(b))/evalf(b,12)}

\subsubsection{Erreur absolue, relative, arrondi 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'{\bf erreur absolue}\index{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\`u la notion
d'erreur {\em relative}~:
\begin{defn}
L'{\bf erreur relative}\index{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}+(2.0^{-53}+1.0) \rightarrow 1 \]
Dans Xcas, il n'y a que 48 bits de mantisse~:
\giacinputbigmath{(2.^(-48)+(2.^(-48)+1.))-1.;  ((2.^(-48)+2.^(-48))+1.)-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$ car le terme
d'erreur est en $O(h^2)$ et non en $O(h)$. Attention
toutefois \`a ne pas prendre $h$ trop petit, sinon $x+h=x$
en flottants et m\^eme si $x+h \neq x$, l'erreur absolue sur
$f(x+h)-f(x-h)$ est (au moins) d'ordre $\varepsilon |f(x)|$, donc l'erreur
relative est d'ordre $\varepsilon/h |f(x)|$. 
Par exemple pour {\tt h=1e-8} le reste est en $O(h^2)$ donc
de l'ordre des erreurs d'arrondi mais l'erreur relative sur
$f(x+h)-f(x-h)$ est d'ordre $\epsilon/h$ largement sup\'erieure
(en flottants double-pr\'ecision). On choisira plut\^ot $h$ tel que
$\epsilon/h$ soit proche de $h^2$, donc de l'ordre de {\tt 1e-5}, qui
fournira une valeur approch\'ee avec une erreur relative de l'ordre
de {\tt 1e-10}.
Exemple~: calcul de la d\'eriv\'ee num\'erique de $\exp(\sin(x))$ en $x=1$
\giacinputbigmath[style="width:800px;height:40px;font-style:large"]{f(x):=exp(sin(x)); seq(taux_accroissement(f(x),1.0-10^(-k),1.0+10^(-k))-f'(1.0),k,1,11)}

Remarquons néanmoins que les erreurs calculées ici sont des majorations
des erreurs réelles (ou si on préfère l'erreur obtenue dans le pire
des cas), statistiquement les erreurs sur les résultats sont moindres,
par exemple si on effectue $n$ calculs susceptibles de provoquer
des erreurs ind\'ependantes suivant une m\^eme loi d'esp\'erance nulle, la moyenne des
erreurs divis\'ee par l'\'ecart-type de la loi
tend vers une loi normale centr\'ee r\'eduite. De mani\`ere plus
d\'eterministe, on a l'in\'egalit\'e de Bienaym\'e-Tchebyshev
\[ P(|X|>\alpha) \leq \frac{n\sigma^2}{\alpha^2}\]
o\`u $X$ est la variable al\'eatoire somme des $n$ erreurs, $\alpha$
l'erreur et $n\sigma^2$ la variance de la somme $n$ erreurs
suppos\'ees ind\'ependantes, cette probabilit\'e tend vers 0 pour $n$
grand si $\alpha$ est d'ordre $n$, et ne tend
pas vers 0 si $\alpha$ est de l'ordre de $\sqrt{n}$.
Exemple~: somme de $n=400$ nombres r\'epartis sur $[-1,1]$ selon la loi
uniforme (repr\'esentant des erreurs), on divise par $\sqrt{n}=$20, 
on effectue plusieurs tirages (par exemple 500) on trace l'histogramme et 
on compare avec la loi normale de moyenne
nulle (l'esp\'erance de la somme) et d'\'ecart-type celui de la loi
uniforme.
\giacinput[style="width:800px;height:40px;font-style:large"]{m:=ranm(400,500,-1..1):;gl_x=-2..2;histogram(sum(m)/20/stddev(uniform,-1,1),-1,0.1); plot(normald(0,1),-2..2)}

Attention, si on effectue la somme de $n$ r\'eels $\sum_j x_j$, 
les erreurs d'arrondis
ne satisfont pas \`a ces hypoth\`eses. En effet, l'erreur d'arrondi
\`a chaque op\'eration est une erreur relative, l'erreur absolue 
correspondante est $\epsilon |x_1+x_2|$ puis $\epsilon |x_1+x_2+x_3|$
puis ... $\epsilon |x_1+x_2+...+x_n|$, que l'on peut majorer par
$$ \epsilon ((n-1)|x_1|+(n-2)|x_2|+...+|x_n||)$$
La majoration de l'erreur d'arrondi d\'epend donc de l'ordre des termes,
on a int\'er\^et \`a sommer en commen\c{c}ant par les termes les
plus petits en valeur absolue.
Voici des programmes C++ illustrant cela
\begin{itemize}
\item
\ahref{http://www-fourier.univ-grenoble-alpes.fr/\home{parisse}/mat406/erreur2.cc}
{uniquement avec des flottants}.
\`A compiler avec la commande \verb|c++ erreur2.cc -o erreur2|, ne pas
utiliser d'option d'optimisation sinon les flottants interm\'ediaires seront
probablement calcul\'e en pr\'ecision plus grande (64 bits de mantisse)
\item
\ahref{http://www-fourier.univ-grenoble-alpes.fr/\home{parisse}/mat406/erreur.cc}{
en comparant avec des flottants multi-pr\'ecision} (n\'ecessite
d'avoir Giac install\'e sous Unix/Linux, se compile avec 
\verb|c++ erreur.cc -lgiac -lgmp -o erreur|).
\end{itemize}
Mais on peut faire mieux, 
il est possible de corriger les erreurs d'arrondi dans une somme
avec le programme suivant pour une liste (on peut bien sur
adapter \`a la somme d'une expression d\'ependant d'une variable
enti\`ere sans stocker de liste)~:
\begin{giaconload}
Somme(l):={
  local x,s,c;
  s:=0.0;
  c:=0.0;
  pour x in l faire
    c += (x-((s+x)-s));
    s += x;
  fpour;
  print(c);
  return s+c;
}:;
\end{giaconload}
En effet, $c$ devrait valoir 0 sans erreurs d'arrondis,
avec les erreurs d'arrondis, on a le premier calcul $s+x$ qui donnera
une erreur oppos\'ee \`a celui du calcul de $s$ \`a la ligne
suivante, le 2i\`eme calcul effectu\'e $(s+x)-s$ donne une erreur
absolue en $\epsilon |x|$ au pire (car c'est une erreur relative
par rapport \`a $(s+x)-s$), 
et la 3i\`eme erreur d'arrondi est n\'egligeable
(puisque la somme vaut 0). On a donc une erreur absolue sur $s+c$
qui est au premier ordre au pire en $O(\epsilon \sum|x_i|)$, 
bien meilleure que
la majoration $\epsilon ((n-1)|x_1|+(n-2)|x_2|+...+|x_n||)$
calcul\'ee pr\'ec\'edemment.


Par exemple
\giacinput{Digits:=14; n:=25000; l:=seq(1/j,j,1,n):; Somme(l);}
\`a comparer avec
\giacinput{S:=sum(1/j,j,1,n):; evalf(S); sum(1. /j,j,1,n);}
(le calcul de $S$ est fait en exact, celui de {\tt sum(1. /j,j,1,n)}
est approch\'e sans correction).

En conclusion, il est souvent tr\`es difficile de calculer une majoration
rigoureuse de l'erreur pour des calculs (sauf les plus simples), et cette
majoration est en g\'en\'eral bien trop pessimiste.
Lorsqu'on doute de la précision d'un calcul, un test peu couteux consiste
à refaire ce calcul en utilisant des flottants en précision plus
grande et tester si le résultat varie en fonction du nombre de chiffres
significatifs utilisés, ou faire varier l\'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 (section suivante)


\subsection{L'arithm\'etique d'intervalle.}\index{intervalle, arithm\'etique}

Certains syst\`emes de calcul formel peuvent manipuler directement
des intervalles r\'eels, par exemple par l'interm\'ediaire de la
bibliothèque C MPFI. Les op\'erations arithm\'etiques sur des
intervalles renvoient alors le meilleur intervalle possible contenant
toutes les valeurs possibles lorsque les op\'erandes parcourent
leurs intervalles respectifs.
Exemple en Xcas (version 1.1.1 et ult\'erieures)~:
\verb|[-1..2]*[-1..2]| renvoie \verb|[-2..4]|.
Attention ici on parcourt toutes les valeurs possibles de 
$xy, \ x \in [-1,2], y \in [-1,2]$. Ce qui est diff\'erent du carr\'e
d'un intervalle ou plus g\'en\'eralement de l'\'evaluation
d'un polyn\^ome en un intervalle, \verb|horner(x^2,[-1..2])|
renvoie ainsi \verb|[0..4]|.

Les fonctions disponibles sont souvent moins riches qu'en
arithm\'etique flottante, le calcul d'une fonction non monotone
sur un intervalle peut s'av\'erer d\'elicat, alors que si la fonction
est monotone, il suffit de calculer l'image des deux bornes
de l'intervalle. Pour les polyn\^omes, Xcas d\'ecompose
les coefficients en deux parties $P=P_+-P_-$ en fonction
du signe, puis utilise la monotonie de $P_+$ et $P_-$ sur
$\mathbb{R}^+$ et $\mathbb{R}^-$ respectivement.

L'arithm\'etique d'intervalle dans $\mathbb{C}$ est beaucoup plus difficile
\`a mettre en oeuvre puisqu'il n'y a plus d'ordre ni de monotonie,
on doit alors s'en remettre \`a des estimations sur les parties
r\'eelles et imaginaires qui ne tiendront pas compte du ph\'enom\`ene
ci-dessus sur la diff\'erence entre $xy, \ x \in [-1,2], y \in [-1,2]$
et $x^2, \ x \in [-1,2]$.




\subsection{Types compos\'es.}
\subsubsection{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\'ecroit \`a chaque it\'eration. Arrive un
moment o\`u le reste est nul, le PGCD est alors l'entier par lequel
on a divis\'e.
Il existe des variantes de cet algorithme un peu plus efficaces
lorsque les nombres sont repr\'esent\'es 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\'erateur et
du d\'enominateur pour l'\'ecrire sous forme irr\'eductible.

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 souvent g\'er\'es par le microprocesseur (ou son
coprocesseur arithm\'etique).

\subsubsection{Les complexes}
Apr\`es les nombres r\'eels, on passe aux 
{\bf nombres complexes}\index{complexe}~: 
on utilise un couple (partie r\'eelle, imaginaire) de fractions (exacts) 
ou de flottants et les r\`egles habituelles sur les complexes.

\subsubsection{Les polyn\^omes}
Apr\`es les nombres, l'objet le plus utilis\'e dans les syst\`emes
de calcul formel est probablement le {\bf polyn\^ome}\index{polynome}, 
toute simplification
d'une expression se ram\`ene \`a un moment donn\'e \`a mettre sous
forme irr\'eductible une fraction de polyn\^omes. Les principales
repr\'esentations possibles sont~:
\begin{itemize}
\item les polyn\^omes \`a 1 variable, repr\'esentation dense, on stocke
la liste des coefficients du polyn\^ome par ordre croissant ou d\'ecroissant
\item les polyn\^omes \`a 1 variable, repr\'esentation creuse, on stocke
des paires coefficients, degr\'e pour les coefficients non nuls
\item les polyn\^omes \`a plusieurs variables, repr\'esent\'e r\'ecursivement
de mani\`ere dense ou creuse (i.e. $P(x_1,...,x_n)$ vu comme polyn\^ome
en $x_n$ \`a coefficients polyn\^omes d\'ependant des variables $x_1,...,x_{n-1}$),
ce sont des cas particuliers des 2 cas pr\'ec\'edents
\item les polyn\^omes \`a plusieurs variables distribu\'es, on stocke
des mon\^omes, qui sont des paires coefficient, liste d'entiers,
la liste repr\'esentant les exposant des variables dans le mon\^ome.
\item la repr\'esentation symbolique (par exemple $x y^2-5x+y^3$)
beaucoup plus difficile \`a manipuler directement
\end{itemize}
Algorithmes de base sur les polyn\^omes~: 
l'\'evaluation 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).

\subsubsection{Calcul symbolique}
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\'epend d'un contexte d'ex\'ecution
et le remplacement du symbole par sa valeur affect\'ee est appel\'e
\'evaluation).

Les {\bf expressions}\index{expression}, par exemple \verb|sin(x)+2*x^2|,
elles peuvent \^etre repr\'esent\'ees par des arbres. L'\'evaluation
d'une expression consiste \`a remplacer les symboles de l'expression
par leur valeur, puis \`a effectuer les op\'erations en tenant compte
de la substitution. Il est parfois souhaitable de ne pas
effectuer certaines op\'erations de substitution, on emp\^eche
l'\'evaluation, explicitement (\verb|''|) ou implicitement
(par exemple l'affectation n'\'evalue 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.

\subsubsection{Listes, s\'equences, tables}
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\'es pour les {\bf vecteurs, matrices}\index{vecteur}\index{matrice}) 
qui sont essentiellement des listes ou des listes de listes de m\^eme taille. 
Les s\'equences sont des listes d'objets ordonn\'es ``non r\'ecursifs''
(ils ne peuvent contenir des s\'equences), alors que les listes
peuvent contenir des listes, sinon il n'y a pas de diff\'erences.
Dans les logiciels de calcul formel, la plupart du temps
les s\'equences se notent en indiquant les \'el\'ements s\'epar\'es par
des virgules. Les listes s'en distinguent par les d\'elimiteurs \verb|[]|.
Il faut prendre garde au fait qu'en g\'en\'eral affecter par exemple
\verb|l[1]:=3;| \`a une variable libre \verb|l| cr\'ee 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\`ere 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\'ees TD: Exercices sur les nombres en base 2,
%sur les changements de base, sur les erreurs.
%Ecriture en base $z$ avec reste sym\'etrique.
%L'addition n'est pas associative dans les flottants : comment
%sommer des s\'eries de terme g\'en\'eral d\'ecroissant pour avoir la meilleure
%pr\'ecision possible. 
%Changement d'origine pour un polynome avec appel de Taylor.


\section{Suites it\'eratives et applications}
R\'esum\'e:\\
Th\'eor\`eme du point fixe, m\'ethode de Newton, convexit\'e. 
Exemple: calcul de valeur approch\'ee de racines carr\'ees,
R\'esolution d'\'equations.

\subsection{Rappel~: suite r\'ecurrente}
Une suite r\'ecurrente \`a un cran est d\'efinie par une valeur
initiale $u_0$ et une relation de r\'ecurrence que l'on va supposer
ne pas d\'ependre explicitement de $n$, donc de la forme
$$ u_{n+1}=f(u_n)$$
On peut repr\'esenter graphiquement une suite de ce type par un graphe
``en toile d'araign\'ee'', par exemple pour les cinq premiers termes
de $f(x)=\frac{x+2}{x+1}$ et $u_0=3$,
\begin{giaconload}
gl_x=-0.5..2.5;plotseq((x+2)/(x+1),[0,0,2],5)
\end{giaconload}
La suite peut converger ou pas, par exemple avec
$f(x)=\tan(x)-0.1$ et $f(x)=\cos(x)$
\giacinput{plotseq(tan(x)-0.1,[0,-1.3,1.3],7);}
\giacinput{plotseq(cos(x),[0,-2,2],5);}

\subsection{Le point fixe}
Soit $f$ une fonction continue sur un intervalle $I=[a,b]$ de $\mathbb{R}$, et 
\`a valeurs dans $I$ (attention \`a bien choisir $I$ pour que l'image
de $I$ par $f$ reste dans $I$).
On s'int\'eresse \`a 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\'erifier
\[ f(l)=l \]
puisque $f$ est continue. On dit que $l$ est un point fixe de $f$.
Ceci am\`ene \`a l'id\'ee d'utiliser ces suites pour r\'esoudre num\'eriquement
l'\'equation $f(x)=x$. Nous allons donner un th\'eor\`eme 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 $f$ est d\'efinie de $I$ dans $I$ et
\[ \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\`erera seront continument
d\'erivables, donc d'apr\`es le th\'eor\`eme des accroissements finis
\[ f(y)-f(x)=f'(\theta) (y-x), \quad \theta \in [x,y] \]
ainsi pour v\'erifier que $f$ est contractante, on \'etudie la valeur absolue
de $f'$ sur $I$, il suffit de montrer que cette valeur absolue
est strictement inf\'erieure \`a un r\'eel $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 \`a 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\'emonstration~: Tout d'abord si $f$ est contractante, on montre \`a partir
de la d\'efinition de la continuit\'e que $f$ est continue. 
Soit $g(x)=f(x)-x$, alors $g$ est continue, positive en $a$ et n\'egative
en $b$, il existe donc $l\in[a,b]$ tel que $g(l)=0$ (th\'eor\`eme des
valeurs interm\'ediaires). Soit $u_n$ une suite d\'efinie 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\'ecurrence \'evidente~:
\[ |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\'eom\'etrique $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 \`a la preuve de la majoration (\(\ref{eq:u_n_l}\)) qui est importante
en pratique car elle donne un test d'arr\^et de calcul des
termes de la suite r\'ecurrente, on \'ecrit 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\'egalit\'e 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\`u le r\'esultat.

{\bf Exemple 1}~: \\
Cherchons une valeur approch\'ee de $\sqrt{2}$ par cette m\'ethode.
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\'erifie que $f(\sqrt{2})=\sqrt{2})$, puis que $f'=-1/(x+1)^2$
donc $f$ d\'ecroit. On va voir si les hypoth\`eses du th\'eor\`eme du point
fixe s'appliquent sur par exemple $[1,2]$. Comme $f$ est d\'ecroissante
$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\'erer la suite \`a partir par exemple de $u_0=1$ et on va converger
vers $\sqrt{2}$ (en s'en rapprochant \`a chaque cran d'un rapport
inf\'erieur \`a $1/4$).
\giacinputmath{f:=(x+2)/(x+1); f1:=factor(f'); x0:=fMin(f1,x=1..2); f1(x=x0);}
\giacinputmath{u:=1.0; pour j de 1 jusque 10 faire u:=f(x=u); print(j,u); fpour; u;}

{\bf Exemple 2}~: \\
Consid\'erons l'\'equation en $x$
\[ x- e \sin(x) =t, \quad e \in [0,1[ \]
c'est l'\'equation du temps utilis\'ee en astronomie pour trouver la
position (plus pr\'ecis\'ement l'angle avec le grand axe de l'ellipse) 
d'une plan\`ete \`a l'instant $t$ 
sur son orbite elliptique ($e$ \'etant l'excentricit\'e
de l'ellipse).
\begin{giaconload}
gl_x=-3.5..3.5;gl_y:=-2.2..2.2; 
E:=ellipse(-1,1,2,color=red); cercle(0,2); theta:=2; 
line(0,slope=tan(theta)); angle(0,5,5*exp(i*theta),"x");
D:=droite(x=2*cos(theta)); M:=inter(D,E,i,legend="M(t)");
E:=ellipse(-1,1,2,color=red); cercle(0,2); theta:=2; 
line(0,slope=tan(theta)); angle(0,5,5*exp(i*theta),"x");
D:=droite(x=2*cos(theta)); M:=inter(D,E,i,legend="M(t)",color=red)
\end{giaconload}
Il n'y a pas de formule exacte permettant de calculer $x$ en fonction de $t$.
Si on a une valeur num\'erique pour $t$, on peut trouver une valeur
num\'erique approch\'ee de $x$ par la m\'ethode du point fixe, en r\'e\'ecrivant
l'\'equation sous la forme 
\[ f(x)=t+e\sin(x) = x \]
On observe que $f$ envoie $\mathbb{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\'eor\`eme s'applique, il suffit de
prendre une valeur initiale dans $[t-e,t+e]$ et d'it\'erer la suite
jusqu'\`a obtenir la pr\'ecision d\'esir\'ee. Par exemple si on veut une
valeur approch\'ee de $x$ \`a $10^{-6}$ pr\`es, il suffira que la diff\'erence
entre deux termes successifs de la suite $u_n$ v\'erifie
\[ |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} \]
Par exemple, pour $t=\pi/2 $ et $E=0.0167$ (la Terre), au bout de 3 it\'erations\\
\giacinput{E:=0.0167; t:=pi/2.; f:=t+E*sin(x); u0:=t; u1:=f(x=u0); u2:=f(x=u1); u3:=f(x=u2); abs(u3-u2)}

Cette m\'ethode n'est pas toujours optimale, car la vitesse de convergence
vers la limite $l$ est dite ``lin\'eaire'', c'est-\`a-dire 
que le temps de calcul pour
avoir $n$ d\'ecimales est proportionnel \`a $n$ (ou
encore il faut effectuer un nombre d'it\'erations
proportionnel \`a $n$, chaque it\'eration faisant
gagner en pr\'ecision de l'ordre du rapport $k$ de contractance). 
En effet, supposons que $f'$ est continue en $l$ et que $0<L=|f'(l)|<1 $.
Il existe alors un intervalle $I=[l-\eta,l+\eta]$ tel que
\[  x \in I \Rightarrow \frac{L}{2} \leq |f'(x)| \leq \frac{1+L}{2} \]
Le th\'eor\`eme des accroissements finis donne alors
\[ |u_{n+1} - l | = |f(u_n)-f(l)| = |f'(\theta)| |u_n-l|,
\quad \theta \in [u_n,l] \]
Si $u_0 \in I$, alors $\theta \in I$ donc $|u_1-l| \leq |u_0-l|$ et
$u_1 \in I$, par r\'ecurrence on a pour tout $n$, $u_n \in I$ 
\[ \frac{L}{2} |u_n-l| \leq |u_{n+1} - l| \leq \frac{1+L}{2} |u_n-l| \]
on a donc par r\'ecurrence 
\[ \left(\frac{L}{2}\right)^n|u_0-l| \leq |u_n-l| \leq  
\left( \frac{1+L}{2} \right)^n|u_0-l|
\]
Donc pour avoir $|u_n-l| \leq \epsilon$ il suffit que 
\[ \left( \frac{1+L}{2} \right)^n|u_0-l| \leq \epsilon \Rightarrow
n \geq \frac{\ln(\frac{\epsilon}{|u_0-l|})}{\ln( \frac{1+L}{2}) } \]
et il faut que
\[ \left(\frac{L}{2}\right)^n |u_0-l| \leq \epsilon
\Rightarrow
n \geq \frac{\ln(\frac{\epsilon}{|u_0-l|})}{\ln( \frac{L}{2}) }
 \]

Si $f$ est suffisamment r\'eguli\`ere,
il existe une m\'ethode plus rapide lorsqu'on est proche de la racine ou lorsque
la fonction a des propri\'et\'es de convexit\'e, c'est la m\'ethode de Newton.
Et m\^eme si Newton n'est pas applicable, une simple dichotomie
peut \^etre plus efficace si la constante de contractance est 
sup\'erieure \`a $1/2$ (y compris pr\'es de la solution de $f(x)=x$).
\begin{giaconload}
def dicho(f,a,b,eps):
    # local c,niter;
    if f(a)*f(b)>=0: 
        return "erreur: f(a)*f(b)>=0"
    while b-a>eps:
        c=(a+b)/2.0
        if f(a)*f(c)>0:
            a=c
        else:
            b=c
    return c
\end{giaconload}
\giacinput{f(x):=cos(x)-x; dicho(f,0.0,1.0,1e-8);}

Cette m\'ethode a toutefois l'avantage de se g\'en\'eraliser en
dimension sup\'erieure \`a 1, contrairement \`a la dichotomie.

\subsection{La m\'ethode de Newton.}
La m\'ethode de Newton est une m\'ethode de r\'esolution de l'\'equation
$f(x)=0$, attention \`a la diff\'erence avec le th\'eor\`eme du point fixe
qui permet de r\'esoudre num\'eriquement $f(x)=x$.
Si $x_0$ est proche de la racine $r$
on peut faire un d\'eveloppement de Taylor \`a l'ordre 1 de la
fonction $f$ en $x_0$~:
\[ f(x)=f(x_0)+(x-x_0)f'(x_0)+O((x-x_0)^2) \]
Pour trouver une valeur approch\'ee de $r$, on ne garde que la partie
lin\'eaire du d\'eveloppement, on r\'esout~:
\[ f(r)=0 \approx f(x_0) + (r-x_0) f'(x_0) \]
donc (si $f'(x_0)\neq 0$)~:
\[ r \approx x_0 -\frac{f(x_0)}{f'(x_0)}\]
Graphiquement, cela revient \`a tracer la tangente \`a la courbe repr\'esentative
de $f$ et \`a chercher o\`u elle coupe l'axe des $x$.
\begin{giaconload}
gl_x=-3..3;
f:=x^2-2; 
G:=plot(f,x=-3..3);
x0:=1.9; M:=element(G,x0);
T:=tangent(G,x0,legend="",color=red);
N:=inter_unique(T,droite(y=0));
legend(-3,"xN="+round(abscisse(N),3))
\end{giaconload}
On consid\`ere donc la suite r\'ecurrente d\'efinie par une valeur $u_0$
proche de la racine et par la relation~:
\[ u_{n+1} = u_n -\frac{f(u_n)}{f'(u_n)}\]
{\bf Remarque~:} C'est un cas particulier de point fixe pour
l'\'equation $g(x)=x$ avec
$$ g(x)=x-\frac{f(x)}{f'(x)}, \quad g'(x)=-\frac{f(x) f'{'}(x)}{f'(x)^2}$$
et donc $g'$ s'annule en un point o\`u $f$ est nulle. On esp\`ere
ainsi avoir une convergence tr\`es rapide pr\`es d'une racine.

Il y a deux th\'eor\`emes importants, l'un d'eux prouve que
si $u_0$ est ``assez proche'' de $r$ alors la suite $u_n$ converge vers $r$,
malheureusement il est difficile de savoir en pratique si on est 
``assez proche'' de $u_0$ pour que ce th\'eor\`eme s'applique. Le second
th\'eor\`eme donne un crit\`ere pratique facile \`a v\'erifier qui assure
la convergence, il utilise les propri\'et\'es de convexit\'e de la fonction.

\begin{theorem} \index{Newton}
Soit $f$ une fonction de classe $C^2$ (2 fois continument d\'erivable)
sur un intervalle ferm\'e $I$. Soit $r$ une racine simple de $f$
situ\'ee \`a l'int\'erieur de $I$
(telle que $f(r)=0$ et $f'(r)\neq 0$). Alors il existe $\varepsilon>0$
tel que la suite d\'efinie 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\'emonstration}~: 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\'eveloppement de Taylor de $f$ en $u_n$ \`a l'ordre 2,
on obtient pour un r\'eel $\theta$
situ\'e 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'{'}|<M$ et $|1/f'|<m$
sur $[r-\varepsilon,r+\varepsilon]$ (c'est toujours possible car
$f'{'}$ et $1/f'$ sont continues au voisinage de $r$ puisque $f'(r)\neq 0$).
Si $u_n$ est dans cet intervalle, alors $\theta$ aussi donc
\[ |u_{n+1}-r| \leq |u_n-r|^2 \frac{Mm}{2} \leq  \frac{|u_n-r|Mm}{2}
|u_n-r|,  \]
On a $|u_n-r| \leq \varepsilon$, on diminue si n\'ecessaire
$\varepsilon$ pour avoir $\varepsilon < 2/(Mm)$, on a alors~:
\[ |u_{n+1}-r| \leq k |u_n-r|, \quad k=\frac{\varepsilon Mm}{2}<1  \]
donc d'une part $u_{n+1}$ est encore dans l'intervalle 
$[r-\varepsilon,r+\varepsilon]$ 
ce qui permettra de refaire le m\^eme raisonnement au rang
suivant, et d'autre part
on a une convergence au moins g\'eom\'etrique vers $r$.
En fait la convergence est bien meilleure
lorsqu'on est proche de $r$ grace au carr\'e dans $|u_n-r|^2$,
plus pr\'ecis\'ement, on montre par r\'ecurrence que
\[ |u_n-r| \leq |u_0 - r|^{2^n} \left( \frac{Mm}{2} \right)^{2^n-1}
\]
il faut donc un nombre d'it\'erations proportionnel \`a $\ln(n)$
pour atteindre une pr\'ecision donn\'ee.

{\bf Remarque~:} ce th\'eor\`eme se g\'en\'eralise sur $\mathbb{C}$ et m\^eme sur $\mathbb{R}^n$.

{\bf Exemple~:} pour calculer $\sqrt{2}$, on \'ecrit l'\'equation $x^2-2=0$
qui a $\sqrt{2}$ comme racine simple sur $I=[1/2,2]$, 
on obtient la suite r\'ecurrente
\[ u_{n+1} = u_n - \frac{u_n^2-2}{2u_n} \]
Si on prend $\eta=1/2$, on a $f'=2x$ et $f'{'}=2$
donc on peut prendre $M=2$ et $m=1$ car $|1/f'|\leq 1$ sur 
$[\sqrt{2}-1/2,\sqrt{2}+1/2]$. On a $2/(mM)=1$, on peut donc 
prendre $\varepsilon=1/2$, la suite convergera pour tout 
$ u_0 \in [\sqrt{2}-1/2,\sqrt{2}+1/2]$.

Plus g\'en\'erallement, on peut calculer une racine $k$-i\`eme d'un r\'eel $a$
en r\'esolvant $f(x)=x^k-a$ par la m\'ethode de Newton.

L'inconv\'enient de ce th\'eor\`eme
est qu'il est difficile de savoir si la valeur de d\'epart qu'on
a choisie se trouve suffisamment pr\`es d'une racine pour que
la suite converge. Pour illustrer le ph\'enom\`ene, 
on peut par exemple colorer les points du plan
complexe en $n+1$ couleurs selon que la suite d\'efinie par la m\'ethode
de Newton converge vers l'une des $n$ racines d'un polyn\^ome de degr\'e
$n$ fix\'e au bout de par exemple 50 it\'erations (la $n+1$-i\`eme couleur
servant aux origines de suite qui ne semblent pas converger).
\ahref{https://xcas.univ-grenoble-alpes.fr/forum/viewtopic.php?f=34&t=2765}{Session
Xcas}

Passons maintenant \`a un crit\`ere tr\`es utile en pratique~:
\begin{defn} (convexit\'e)\index{convexe}\\
Une fonction $f$ continument d\'erivable sur un intervalle $I$ de $\mathbb{R}$
est dite convexe si son graphe est au-dessus de la tangente en tout point
de $I$.
\end{defn}
Il existe un crit\`ere simple permettant de savoir si
une fonction de classe $C^2$ est convexe~:
\begin{theorem}
Si $f$ est $C^2$ et $f'{'} \geq 0$ sur $I$ alors $f$ est convexe.
\end{theorem}
{\bf D\'emonstration}~:\\
L'\'equation de la tangente au graphe en $x_0$ est 
\[ y=f(x_0)+f'(x_0)(x-x_0) \]
Soit
\[ g(x)=f(x)-(f(x_0)+f'(x_0)(x-x_0)) \]
on a~:
\[ g(x_0)=0, \quad g'(x)=f'(x)-f'(x_0), \quad g'(x_0)=0, 
\quad g'{'}=f'{'} \geq 0 \]
donc $g'$ est croissante, comme $g'(x_0)=0$, $g'$ est n\'egative
pour $x<x_0$ et positive pour $x>x_0$, donc $g$ est d\'ecroissante
pour $x<x_0$ et croissante pour $x>x_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\`eme th\'eor\`eme sur la m\'ethode 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\'ethode de Newton
\[ u_{n+1} = u_n -\frac{f(u_n)}{f'(u_n)},  \]
est d\'efinie, d\'ecroissante, minor\'ee par $r$ et converge vers 
$r$. De plus
\[ 0 \leq u_n -r \leq \frac{f(u_n)}{f'(r)} \]
\end{theorem}

{\bf D\'emonstration}~:\\
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\'eduit que $f>0$ sur $]r,b]$ donc $u_{n+1} \leq u_n$.
Comme la courbe repr\'esentative 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\'ecroissante minor\'ee 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\`ere $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$.

{\bf Remarque~:}
Lorsque $u_n$ est proche d'une racine, le calcul approch\'e de $f(u_n)$
peut poser des probl\`emes \`a cause de compensations entre deux
valeurs proches que l'on soustrait. Par exemple pour $f(x)=x^2-2$,
si $x$ est proche de $\sqrt{2}$, $x^2$ sera proche de 2, et le calcul
de la diff\'erence  va perdre \'enorm\'ement de pr\'ecision relative.
\giacinputbigmath{
f(x):=x^2-2;
g:=id-f/f';
u:=2.0; l:=[[u,f(u)]];
pour j de 1 jusque 6 faire u=g(u); l.append([u,f(u)]); fpour:;
l;
}
On peut contourner ce probl\`eme en utilisant des flottants
multip\'recision, ci-dessus en posant par exemple \verb|u=evalf(u,40)|.
Dans le cas de la recherche de
racines de polyn\^omes, on peut utiliser des rationnels (\verb|u=2|) mais un autre probl\`eme
se produit~: la taille du num\'erateur et du d\'enominateur va en gros
doubler \`a chaque it\'eration! Pour \'eviter cela, on peut remplacer
par un rationnel proche de taille de num\'erateur/d\'enominateur plus
petite par arrondi \`a une puissance de 2 d\'ependant de la
pr\'ecision attendue de $u_n$.
%% Id\'ees TD: diff\'erentes suites it\'eratives approchant $\sqrt{3}$,
%% par exemple trouver (a*x+b)/(c*x+d) ayant sqrt(3) comme point fixe
%% vitesse de convergence th\'eorique.
%% Recherche de suite it\'erative pour les racines $n$-i\`eme.

%% Id\'ees TP: observation pratique de la vitesse de convergence,
%% en pr\'ecision double ou illimit\'ee. Exemple de cas o\`u Newton marche
%% et ne marche pas, utilisation du pr\'efacteur.

\section{D\'eveloppement de Taylor, s\'eries enti\`eres, fonctions usuelles}
R\'esum\'e: S\'eries enti\`eres. Calcul des fonctions
transcendantes usuelles. 

Soit $f$ une fonction ind\'efiniment d\'erivable sur un intervalle $I$ de $\mathbb{R}$
et $x_0 \in I$. On peut alors effectuer le d\'eveloppement de Taylor de $f$
en $x_0$ \`a 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 \'egale \`a $f(x)$ et si on peut facilement
majorer la diff\'erence entre $f(x)$ et $T_n(f)(x)$. Si c'est le
cas, on pourra utiliser $T_n(f)(x)$ comme valeur approch\'ee de $f(x)$.

On peut parfois r\'epondre \`a ces questions simultan\'ement en regardant le
d\'eveloppement de Taylor de $f$ avec reste~
\begin{thm}
Si $f$ est de classe $C^{n+1}$ sur $[x_0,x]$ (ou $[x,x_0]$),
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)!} \]
\end{thm}
La preuve se fait en consid\'erant la fonction
$$g(t) := T_{n,x}(f)(t)+C(x-t)^{n+1}-f(t), \quad
T_{n,x}(f)(t)= f(x) + (t-x) f'(x) + ... + (t-x)^n \frac{f^{[n]}(x)}{n!}$$
qui v\'erifie $g(x)=0$ ainsi que toutes ses dérivées jusqu'à l'ordre $n$.
On choisit $C$ pour que $g(x_0)=0$ et on
déduit que la d\'eriv\'ee s'annule en $\theta_0$ 
entre $x_0$ et $x$, puis la dérivée seconde s'annule
en $\theta_1$ entre $\theta_0$ et $x$ donc entre $x_0$ et $x$,
etc. ce qui donne la valeur de $C$ pour la dérivée $n+1$-ième de $g$.

C'est le cas pour la fonction exponentielle que nous allons
d\'etailler, ainsi que les fonctions sinus et cosinus.

\subsection{La fonction exponentielle}\index{exp}
Soit $f(x)=\exp(x)$ et $x_0=0$, la d\'eriv\'ee $n$-i\`eme 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\'egatif,
$|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\'eel
\[ 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!} \]
\begin{giaconload}
gl_x=-1.5..1.5,gl_y=-4..4
,plot([exp(x),1+x,1+x+x^2/2,1+x+x^2/2+x^3/6,
1+x+x^2/2+x^3/6+x^4/24],x=-1.5..1.5,
color=[black,blue,cyan,green,red])
\end{giaconload}

Comment en d\'eduire une valeur approch\'ee 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\'erieur \`a 
l'erreur absolue souhait\'ee, le plus t\^ot \'etant le mieux pour des
raisons d'efficacit\'e et pour \'eviter l'accumulation d'erreurs
d'arrondi. 
Si on veut connaitre $e^x$ \`a une erreur relative $\varepsilon$ donn\'ee
(par exemple $\varepsilon=2^{-53}$ pour stocker le r\'esultat 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\^eter la sommation
lorsque le terme suivant est plus petit que $\varepsilon$.

On observe que plus $x$ est grand, plus $n$ devra
\^etre grand pour r\'ealiser le test d'arr\^et, ce qui est facheux
pour le temps de calcul.
De plus, le r\'esultat final peut \^etre petit alors que les termes
interm\'ediaires calcul\'es dans la somme peuvent \^etre grands, ce qui
provoque une perte de pr\'ecision relative, par exemple si on
veut calculer $e^{-10}$ ou plus g\'en\'erallement l'exponentielle
d'un nombre n\'egatif de grande valeur absolue.

Exercice~: combien de termes faut-il calculer dans le d\'eveloppement
de l'exponentielle de -10 pour que le reste soit plus petit
que $2^{-53}$~? Quel est la valeur du plus grand terme rencontr\'e dans
la suite~? Quelle est la perte de pr\'ecision relative occasionn\'e
par cette m\'ethode de calcul~?

On peut utiliser les propri\'et\'es de la fonction exponentielle
pour \'eviter ce probl\`eme. Pour les nombres n\'egatifs, on peut
utiliser l'\'equation $e^{-x}=1/e^x$ (ne change pas l'erreur relative). 
Pour les grands r\'eels, on peut utiliser $e^{2x}=(e^x)^2$
(multiplie par 2 l'erreur relative).
On peut aussi, si on connait une valeur approch\'ee
de $\ln(2)$, effectuer la division euclidienne de $x$ par $\ln(2)$ 
avec reste sym\'etrique~:
\[ x = a \ln(2) + r, \quad a \in \mathbb{Z}, |r| \leq \frac{\ln(2)}{2} \]
puis si $r$ est positif, on somme la s\'erie de $T(f)(r)$, si $r$
est n\'egatif, 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'\'etant pas connu exactement,
on commet une erreur d'arrondi absolu sur $r$ d'ordre $a \eta$,
o\`u $\eta$ est l'erreur relative sur $\ln(2)$,
il faut donc ajouter une erreur d'arrondi relative de $x \eta /\ln(2)$
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\'ethode de Horner pour \'evaluer
$T_n(f)(r)$ (car elle commence par sommer les termes de plus haut degr\'e
qui sont justement les plus petits termes de la somme). 
Les coprocesseurs arithm\'etiques qui impl\'ementent la fonction exponentielle
ont un format de repr\'esentation 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\'ecision,
le calcul de $e^x$ par cette m\'ethode entraine donc seulement une erreur
relative d'arrondi au plus proche sur le r\'esultat converti en double 
(donc de $2^{-53}$).

Notons que en g\'en\'eral $x$ lui-m\^eme
a d\'ej\`a \'et\'e arrondi ou n'est connu qu'avec une pr\'ecision 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\'erer mieux qu'une erreur relative de 
$e^{\varepsilon |x|}-1$ sur l'exponentielle de $x$. Si $\varepsilon x$ est petit
cette erreur relative (impossible \`a \'eviter, quel que soit
l'algorithme utilis\'e 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\'ee peut \^etre tr\`es \'eloign\'ee de la valeur
r\'eelle! Notons que pour les double, il y aura dans ce cas d\'ebordement 
soit vers l'infini soit vers 0
(par exemple si $x$ est sup\'erieur \`a 709, l'exponentielle renvoie infini).

Exercice~: refaire les m\^emes 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\'ethode a toutefois ces limites, car il peut devenir impraticable
de calculer la d\'eriv\'ee $n$-i\`eme d'une fonction (par exemple avec $\tan(x)$),
et encore plus de la majorer. D'o\`u l'int\'er\^et de d\'evelopper une th\'eorie
des fonctions qui sont \'egales \`a leur d\'eveloppement de Taylor \`a l'infini
d'une part, et d'avoir d'autres m\'ethodes pour majorer le reste, nous
pr\'esentons ici le cas des s\'eries altern\'ees.

Remarque~: pour calculer la fonction exponentielle en optimisant le nombre
d'op\'erations \`a effectuer, on utilise des approximants de Pad\'e,
ce sont des fractions rationnelles qui ont le d\'eveloppement de
Taylor souhait\'e en $x=0$, et pr\'esentent des sym\'etries permettant
de les calculer efficacement.
\giacinputmath{P:=pade(exp(x),x,10,6); series(P,x=0,10);}
Le programme suivant permet de calculer l'exponentielle sur
$[-\ln(2)/2,\ln(2)/2]$ en 12 op\'erations
avec une pr\'ecision proche de \verb|1e-15|
\begin{giaconload}
def f(x): 
    X=x*x
    q=((X+420)*X+15120)*x
    r=(30X+3360)*X+30240
    return (r+q)/(r-q)
\end{giaconload}
\giacinput{f(0.34); exp(0.34)}
En ajoutant une division par $\ln(2)$ avec arrondi  \`a l'entier le
plus proche, on a l'exponentielle en 16 op\'erations arithm\'etiques
de base.

\subsection{S\'eries enti\`eres.}\index{serie entiere}
Les s\'eries de type prendre la limite lorsque $n$ tend vers
l'infini du d\'eveloppement 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\'eresser plus g\'en\'erallement \`a $\sum_{n=0}^\infty a_n x^n$
lorsque $a_n$ est un complexe quelconque, 
c'est ce qu'on appelle une s\'erie enti\`ere, on peut aussi les voir comme
des polyn\^omes g\'en\'eralis\'es. 

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}} {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}} {1-|\frac{x}{x_0}|} 
\end{equation}
\end{theorem}

On en d\'eduit qu'il existe un r\'eel
positif $R\geq 0$ \'eventuellement \'egal \`a $+\infty$ tel que 
la s\'erie converge (la limite de la somme jusqu'\`a l'infini existe) 
lorsque $|x|<R$ et n'existe pas lorsque
$|x|>R$, ce r\'eel est appel\'e {\bf rayon de convergence} de la s\'erie.
Par exemple ce rayon vaut $+\infty$ pour l'exponentielle, le sinus
ou le cosinus. Il
est \'egal \`a 1 pour la s\'erie g\'eom\'etrique $\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 deux s\'eries ont un rayon de convergence non nul, alors on
peut effectuer leur somme, leur produit comme des polyn\^omes et la
s\'erie somme/produit a un rayon de convergence au moins \'egal au plus
petit des 2 rayons de convergence des arguments. On peut inverser une s\'erie
enti\`ere non nulle en 0 en appliquant 
\[ (1+x)^{-1} = 1-x+x^2-x^3+... \]
et on obtient une s\'erie enti\`ere de rayon de
convergence non nul. On peut aussi composer deux s\'eries enti\`eres
$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\'eriver
et int\'egrer une s\'erie enti\`ere terme \`a terme dans son rayon de convergence.

On dit qu'une fonction est d\'eveloppable en s\'erie enti\`ere en 0 si
elle est \'egale \`a son d\'eveloppement de Taylor en 0 somm\'e jusqu'en l'infini
dans un disque de centre 0 et de rayon non nul. Les fonctions
exponentielle, sinus, cosinus sont donc d\'eveloppables en s\'erie enti\`ere en 0.
La fonction tangente \'egalement car le d\'enominateur 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\'eveloppable en s\'eries enti\`eres
pour tout $\alpha \in \mathbb{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\'erie 
g\'eom\'etrique 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\'egrant par rapport \`a $x$, on obtient que $\ln(1+x)$ est d\'eveloppable
en s\'erie enti\`ere 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\`ere analogue le d\'eveloppement en s\'erie enti\`ere
de $\arctan(x)$ en iint\'egrant celui de $1/(1+x^2)$, de m\^eme pour $\arccos(x)$
et $\arcsin(x)$ en int\'egrant 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\'epondre
\`a la question o\`u arr\^ete-t-on la somme pour obtenir une pr\'ecision donn\'ee? 
Dans le cas de $\ln(1+x)$,
on pourrait r\'epondre comme avec l'exponentielle en majorant la d\'eriv\'ee
$n+1$-i\`eme, mais ce n'est plus faisable pour $\arctan, \arcsin, \arccos$.
On va donner un autre crit\`ere qui ne n\'ecessite pas
de calculer cette d\'eriv\'ee mais utilise l'alternance des signes
dans la somme.

\subsection{S\'erie altern\'ee}\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.
\item En fait pour calculer le logarithme \`a la pr\'ecision usuelle
pour $x\in[1/2,2]$, il existe une autre optimisation,
r\'esoudre $x=\frac{1+t}{1-t}$ et utiliser
le d\'eveloppement en s\'eries 
$$ \ln(\frac{1+t}{1-t})=2\sum_{k=0}^\infty \frac{t^{2k+1}}{2k+1}$$
pour $|t| \leq 1/3$.
Cette m\'ethode n\'ecessite d'aller jusque $k=14$ pour majorer le
reste en \verb|1e-16|, ce qui donne une m\'ethode de calcul du log
en une trentaine d'op\'erations.
\end{itemize}

Nous sommes donc en mesure de calculer pr\'ecis\'ement le logarithme
$\ln(1+x)$ pour $|x|<1/2$. Pour calculer $ln$ sur $\mathbb{R}^+$,
on se ram\`ene \`a $[1/2,1]$ en utilisant l'\'ecriture
mantisse-exposant.

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~:\\
\giacinput{a:=sum((1/2)^k/k,k=1..80)}\\
puis la division en base 2 avec 81 bits de pr\'ecision
\giacinput{q:=iquo(numer(a)*2^81,denom(a)); evalf(q*2^(-81))}

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\'efinies par des int\'egrales (cas de nombreuses
fonctions sp\'eciales).

\subsubsection{Exemple~: la fonction d'erreur (error fonction, {\tt erf})}
Cette fonction est d\'efinie \`a une constante multiplicative pr\`es par~:
\[ f(x)=\int_0^x e^{-t^2} \ dt \]
On peut d\'evelopper en s\'eries enti\`eres l'int\'egrand (rayon
de convergence $+\infty$), puis int\'egrer terme
\`a terme, on obtient
\[ f(x)= \sum_{n=0}^{+\infty} (-1)^n \frac{x^{2n+1}}{n! (2n+1)}\]
Ce d\'eveloppement converge tr\`es rapidement pour $|x|\leq 1$.
Par contre, pour $|x|$ grand, il faut calculer beaucoup de termes
avant que le reste soit suffisamment petit pour \^etre n\'egligeable,
et certains termes interm\'ediaires sont grands, ce qui provoque
une perte de pr\'ecision qui peut rendre le r\'esultat calcul\'e
compl\`etement faux. Contrairement \`a la fonction exponentielle,
il n'y a pas de possibilit\'e de r\'eduire l'argument \`a une plage
o\`u la s\'erie converge vite. Il faut donc
\begin{itemize}
\item soit utiliser des flottants multipr\'ecision, avec une pr\'ecision
augment\'ee de la quantit\'e n\'ecessaire pour avoir un r\'esultat fiable
\item soit, pour les grandes valeurs de $x$, utiliser un d\'eveloppement
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\'eveloppement asymptotique
s'obtient par exemple en changeant de variable
$u=t^2$ et en effectuant des int\'egrations par parties r\'ep\'et\'ees
en int\'egrant $e^{-u}$ et en d\'erivant $u^{-1/2}$ et ses d\'eriv\'ees
successives. Ce type de d\'eveloppement asymptotique a la propri\'et\'e
inverse du d\'eveloppement en 0: les termes successifs commencent
par d\'ecroitre avant de croitre et de tendre vers l'infini. Il faut
donc arr\^eter le d\'eveloppement \`a un rang donn\'e (d\'ependant de $x$)
et il est impossible d'obtenir une pr\'ecision meilleure pour cette
valeur de $x$ par un d\'eveloppement asymptotique (on parle parfois de
d\'eveloppement des astronomes).
\end{itemize}

{\bf Exercice~}: donner une valeur approch\'ee de $f(1)$ \`a $1e-16$ pr\`es.
Combien de termes faut-il calculer dans la somme pour trouver
une valeur approch\'ee de $f(7)$ \`a $1e-16$ pr\`es~? Comparer la valeur
de $f(7)$ et la valeur absolue du plus grand terme de la s\'erie,
quelle est la perte de pr\'ecision relative si on effectue les
calculs en virgule flottante~? Combien de chiffres significatifs
faut-il utiliser pour assurer une pr\'ecision finale de 16 chiffres
en base 10~? Calculer le d\'eveloppement asymptotique en l'infini
et d\'eterminer un encadrement de $f(7)$ par ce d\'eveloppement. Combien
de termes faut-il calculer pour d\'eterminer $f(10)$ \`a $1e-16$ pr\`es
par le d\'eveloppement asymptotique et par le d\'eveloppement en s\'eries~?
Quelle est la meilleure m\'ethode 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|<C$, alors le reste de la
s\'erie enti\`ere est major\'e par $|a_nx^n|/(1-|Cx|)$ lorsque
$|x|<1/C$, on peut alors calculer des valeurs approch\'ees
de la fonction solution \`a la pr\'ecision souhait\'ee en utilisant
le d\'eveloppement en s\'eries enti\`eres.

\subsubsection{Exemple~: fonctions de Bessel d'ordre entier}
Soit $m$ un entier positif fix\'e, on consid\`ere l'\'equation
diff\'erentielle
\[ x^2 y'{'} + x y' + (x^2-m^2)y=0 \]
dont on cherche une solution s\'erie enti\`ere 
$y=\sum_{k=0}^\infty a_k x^k $. En remplacant dans l'\'equation, si
$x$ est dans le rayon de convergence de la s\'erie (rayon suppos\'e
non nul), on obtient
\[ 
\sum_{k=0}^\infty k(k-1)a_k x^k + \sum_{k=0}^\infty k a_k x^k 
+ \sum_{k=0}^\infty (x^2-m^2) a_k x^k =0
\]
soit encore
\begin{eqnarray*}
0 &=& \sum_{k=0}^\infty (k^2-m^2+x^2) a_k x^k  \\
 &=& -m^2 a_0 + (1-m^2)a_1 x + \sum_{k=2}^\infty [(k^2-m^2) a_k +a_{k-2}]x^k 
\end{eqnarray*}
Par exemple, prenons le cas $m=0$. On a alors $a_0$ quelconque, $a_1$
nul et pour $k\geq 2$
\[ a_k = - \frac{a_{k-2}}{k^2}\]
Donc tous les $a$ d'indice impair sont nuls. Les pairs sont non nuls
si $a_0\neq 0$, et ils sont de signe altern\'e.
Soit $x$ fix\'e, on observe que pour $2k > |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\'eveloppements asymptotiques et s\'eries divergentes}
Un d\'eveloppement asymptotique est une g\'en\'eralisation d'un d\'eveloppement de Taylor, par exemple
lorsque le point de d\'eveloppement est en l'infini. De nombreuses fonctions ayant
une limite en l'infini admettent un d\'eveloppement asymptotique en l'infini, mais ces
d\'eveloppements sont souvent des s\'eries qui semblent commencer par converger
mais sont divergentes. Ce type de d\'eveloppement s'av\`ere n\'eanmoins tr\`es utile lorsqu'on
n'a pas besoin d'une trop grande pr\'ecision sur la valeur de la fonction.

Nous allons illustrer ce type de d\'eveloppement sur un exemple, la fonction 
exponentielle int\'egrale, d\'efinie \`a une constante pr\`es par
\[ f(x)=\int_x^{+\infty} \frac{e^{-t}}{t} \ dt \]
On peut montrer que l'int\'egrale existe bien, car l'int\'egrand est positif et inf\'erieur \`a $e^{-t}$
(qui admet $-e^{-t}$ comme primitive, cette primitive ayant une limite en $+\infty$).
Pour trouver le d\'eveloppement asymptotique de $f$ en $+\infty$, on effectue des int\'egrations
par parties r\'ep\'et\'ees, en int\'egrant l'exponentielle et en d\'erivant 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\`u
\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\'eveloppement en s\'eries est divergent puisque pour $x>0$ fix\'e et $n$ tendant vers l'infini
\[ \lim_{n\rightarrow +\infty} \frac{n!}{x^{n+1}} = +\infty\]
mais si $x$ est grand, au d\'ebut la s\'erie semble converger, de mani\`ere tr\`es rapide~:
\[ \frac{1}{x} >> \frac{1}{x^2} >> \frac{2}{x^3} \]
On peut utiliser $S(x)$ comme valeur approch\'ee 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\'eries altern\'ees, 
l'erreur relative est inf\'erieure
\`a la valeur absolue du dernier terme somm\'e divis\'e par $e^{-x}/x$. 
Pour $x$ fix\'e assez grand, il 
faut donc trouver un rang $n$, s'il en existe un, 
tel que $(n+1)!/x^{n+1}<\epsilon$ o\`u
$\epsilon$ est la pr\'ecision relative que l'on s'est fix\'ee.
Par exemple, si $x\geq 100$, $n=11$ convient pour 
$\epsilon=12!/100^{12}=5e-16$ (\`a peu pr\`es
la pr\'ecision 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\'e 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\'efinit $\gamma$ comme sa limite. Malheureusement, la convergence
est tr\`es lente et cette d\'efinition n'est pas applicable pour obtenir la valeur
de $\gamma$ avec une tr\`es grande pr\'ecision.
Il y a un lien entre $\gamma$ et la fonction exponentielle int\'egrale, plus pr\'ecis\'ement
lorsque $x\rightarrow 0$, $f(x)$ admet une singularit\'e en $-\ln(x)$,
plus pr\'ecis\'ement $f(x)+\ln(x)$
admet un d\'eveloppement en s\'eries (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\'emonstration condens\'ee) 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\'ee de $\gamma$, il suffit donc de prendre un $x$ assez grand
pour pouvoir calculer $f(x)$ par son d\'eveloppement asymptotique \`a la pr\'ecision requise,
puis de calculer l'int\'egrale du membre de droite par le d\'eveloppement en s\'eries en $x=0$
(en utilisant une pr\'ecision interm\'ediaire plus grande puisque ce d\'eveloppement en s\'eries
va sembler diverger au d\'ebut 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\^eme mani\`ere avec le d\'eveloppement 
en s\'eries et asymptotique
de la fonction sinus int\'egral (on remplace exponentielle par sinus dans
la d\'efinition de $f$) et l'\'egalit\'e (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\'egration par parties, cette fois en int\'egrant $1/t$
et en d\'erivant l'exponentielle (moins 1 dans la premi\`ere int\'egrale).
\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\'egrale, on utilise
l'\'egalit\'e (qui se d\'emontre par r\'ecurrence en faisant une 
int\'egration par parties)~:
\[ n!= \int_0^{+\infty}t^n e^{-t} \ dt \]
On va \`a nouveau int\'egrer par parties,
on int\`egre un facteur multiplicatif 1 
et on d\'erive l'int\'egrand, on simplifie, puis
on int\`egre $t$ et on d\'erive 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\`u
\[ I_n=\int_0^{+\infty} \frac{t^n}{n!} e^{-t} (\ln(t)-\ln(n)) \ dt \]
Pour d\'eterminer $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\^eme 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\'egrales
sont \'equivalentes \`a leur valeur sur un petit intervalle autour de $u=1$, point o\`u l'argument
de l'exponentielle est maximal, 
et comme l'int\'egrand du num\'erateur a une amplitude $\ln(u)$ qui s'annule en $u=1$, 
il devient n\'egligeable devant le d\'enominateur. 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\^omes~: 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 $\mathbb{R}$ ou $\mathbb{C}$ ou $\mathbb{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\'es ci-dessus.
Au rang 2, $Q_0$ est le quotient euclidien de $A$ par $B$ (fonction
\verb|quo|) donc $x$, d'o\`u
\[ 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 \]
\giacinputmath{egcd(x^3-1,x^2+1)}

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\'esoudre
\[ (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\'ee (de degr\'es 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)\]\\
\giacinputmath{abcuv(x^3-1,x^2+1,2x^2)}

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\'eriv\'ee du d\'enominateur au num\'erateur
pour \'eliminer 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\'ecomposer
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$. \\
\giacinputmath{pade(exp(x),x,10,6)}\\
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\'ee est la m\^eme, sauf qu'on remplace l'int\'egrale
par la transform\'ee de Laplace inverse (et les formules donnant
la transform\'ee 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 \`a l'examen)
\item le calcul du terme d'ordre $n$ du d\'eveloppement de Taylor
en 0 d'une fraction rationnelle. On d\'ecompose, et on se ram\`ene \`a des
s\'eries dont le terme g\'en\'eral est connu, comme $(a+x)^{-n}$.
Par exemple pour connaitre le d\'eveloppement de $1/(x^2-3x+2)$,
on factorise le d\'enominateur $1/((x-1)(x-2))$, on d\'ecompose
\[ \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\'eveloppe, le terme d'ordre $n$ est donc $1-(1/2)^ {n+1}$.
%\item Le calcul de sommes de s\'eries 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\'egrale
\[ \int \frac{1}{(x-1)(x^2+1)} \]
en utilisant l'identit\'e de B\'ezout pour d\'ecomposer la fraction rationnelle.
Trouver \`a l'aide de cette d\'ecomposition le terme d'ordre
$n$ du d\'eveloppement de Taylor de la fraction \`a int\'egrer, v\'erifier
avec un logiciel de calcul formel que les termes d'ordre 0 \`a 3 sont
corrects.\\
\giacinputmath{partfrac(1/(x-1)/(x^2+1))}

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$.\\
\giacinputmath{egcd(x^2+y^2-9,x^2+2y^2-2x*y-7)}\\
Cette m\'ethode se syst\'ematise, le polynome obtenu par \'elimination
d'une variable est appel\'e r\'esultant. \\
\giacinputmath{resultant(x^2+y^2-9,x^2+2y^2-2x*y-7,x)}\\
%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}
\giacinpumath{P:=x^4-1; gcd(P,P')}

On a aussi le r\'esultat suivant~:
\begin{prop}
Soit $P$ un polyn\^ome \`a coefficients r\'eels ou complexes, et soit
$Q=P/$pgcd$(P,P')$. Alors $P$ et $Q$ ont les m\^emes racines et
les racines de $Q$ sont simples.
\end{prop}
Preuve~: si $P(z)=C\prod_i (z-z_i)^{m_i}$ avec $z_i$ les racines
distinctes de $P$, alors
pgcd$(P,P')=\prod_i (z-z_i)^{m_i-1}$
et $Q=\prod_i (z-z_i)$.\\
\giacinputmath{P:=x^11-4*x^9-x^8+6*x^7+4*x^6-4*x^5-6*x^4+x^3+4*x^2-1; quo(P,gcd(P,P'))}

Ce r\'esultat est tr\`es utile si le polyn\^ome donn\'e est
\'a coefficients exacts, car les m\'ethodes num\'eriques
d'approximation de racines ne fonctionnent bien que pour
des polyn\^omes \`a racines simples.

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{giaconload}
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{giaconload}
\giacinputmath{yun(x^11-4*x^9-x^8+6*x^7+4*x^6-4*x^5-6*x^4+x^3+4*x^2-1)}
L'instruction \verb|sqrfree| ou \'equivalente
de votre logiciel de calcul formel effectue cette d\'ecomposition.
\giacinputmath{sqrfree(x^11-4*x^9-x^8+6*x^7+4*x^6-4*x^5-6*x^4+x^3+4*x^2-1)}


\subsubsection{Factorisation dans $\mathbb{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 $\mathbb{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 $\mathbb{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 $\mathbb{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$.
Pour des polynômes de degrés grands, cela ne suffit pas,
il existe des méthodes plus efficaces, Maehly, Durand-Kerner-Weierstrass,
Aberth, ... Ainsi la méthode d'Aberth améliore la précision
de racines approchées $z_1,...,z_n$ de $P$ en effectuant la méthode
de Newton sur $P(x)/\prod_{j \neq i}(x-z_j)$, qui est proche de $x-z_i$
sauf tout près des racines approchées.

Si $P$ est \`a coefficients exacts, pour \'eviter les erreurs de
calculs en approch\'e, on remplace $r$ par un rationnel proche et on
effectue ces it\'erations de Newton en calcul exact (entre deux
it\'erations, pour des raisons d'efficacit\'e, on peut remplacer $r$
par un rationnel proche \`a la pr\'ecision attendue \`a cette \'etape
mais dont le num\'erateur et d\'enominateur ont une \'ecriture moins longue).

\subsubsection{Localisation d'une racine complexe pr\`es d'une racine approch\'ee}
On peut montrer \`a post\'eriori des estimations sur la
distance entre une racine approch\'ee et la racine la plus proche
d'un polyn\^ome, plus pr\'ecis\'ement cette distance est inf\'erieure
ou \'egale au degr\'e du polyn\^ome multipli\'e par le module
de $P/P'$ en la racine approch\'ee.

\begin{prop}
Le disque de centre $z$ et de rayon degre$(P) |P(z)/P'(z)|$
contient au moins une racine de $P$.
\end{prop}
En effet
\[ \frac{P'}{P}(z)= \sum_{k} \frac{1}{z-z_k}\]
donc si $|z-z_k| > $degre$(P) |P(z)/P'(z)|$ pour toutes les racines $z_k$ alors
$|\frac{P'}{P}(z)| > |P'(z)/P(z)|$ contradiction.

En combinant des m\'ethodes d'approximation num\'erique des racines
d'un polyn\^ome et ce r\'esultat on peut ainsi localiser tr\`es
pr\'ecis\'ement
des racines complexes de polyn\^ome. C'est un cas o\`u une m\'ethode
hybride exacte-approch\'ee donne des r\'esultats excellents!

Remarque~: pour des polyn\^omes de degr\'e $n$ d\'epassant la dizaine,
on utilise des m\'ethodes permettant de trouver toutes les racines \`a
la fois (diagonalisation num\'erique de la matrice companion) plut\^ot
que la m\'ethode de Newton qui les d\'etermine une par une. 
Lorsque le degr\'e d\'epasse quelques
centaines, il devient n\'ecessaire d'utiliser des approximations
multi-pr\'ecision ou des r\'esutats de fonction \`a variable complexes.

\subsubsection{Localisation des racines r\'eelles~: Sturm}
\index{factorisation}
Pour factoriser un polyn\^ome \`a coefficients r\'eels, on commence
par le factoriser dans $\mathbb{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 $\mathbb{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 $\mathbb{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.\\
\giacinputmath{sturmab(x^3+x+1,-2,2);}

{\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 $0<i<k$, alors en $b$ on a soit ...,-,0,+,... soit
...,+,0,-,... . Regardons le premier cas (le deuxi\`eme cas se traite
de mani\`ere analogue), pour $x$ proche de $b$, on va avoir
...,-,-,+,... ou ...,-,+,+,... dans les 2 cas la contribution au 
nombre de changements de signe est constant (\'egal \`a 1).

Comme $A_k$ est constant, seules les racines de $A_0=P$
sont susceptibles de faire varier $s$. Comme $A_1=P'$, le sens de
variations de $A_0$ au voisinage d'une racine de $A_0$ est d\'etermin\'e
par le signe de $A_1$, donc lorsque $x$ augmente en traversant une
racine $r$ de $P$, il y a deux possibilit\'es soit $P$ est croissant
et on passe de -,+,... \`a +,+,..., soit $P$ est d\'ecroissant et
on passe +,-,... \`a -,-,.... Dans les deux cas, on diminue 
$s$ d'une unit\'e.

{\bf Application}~:\\
Si il n'existe pas de racines r\'eelles dans un intervalle donn\'e,
alors le polyn\^ome garde un signe constant sur cet intervalle,
que l'on peut d\'eterminer en calculant la valeur du polyn\^ome
en un point de cet intervalle. On peut ainsi \'etablir
dans certains cas que la m\'ethode de Newton pour trouver une racine
d'un polyn\^ome convergera.

Par exemple pour le polyn\^ome $P=3x^5-10x^3+30x^2-x-45$, on a
$P'{'}=60(x^3-x+1)$, est positif sur $\mathbb{R}^+$ (exercice~: calculer
la suite de Sturm correspondante pour le v\'erifier). On v\'erifie
que $P(1)<0$ et $P'(1)>0$ donc il existe une racine $r>1$ telle
que $P'(r)>0$, toute valeur de d\'epart de Newton sup\'erieure
\`a $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{Localisation des racines r\'eelles~: r\`egle des signes de Descartes\index{Descartes}\index{regle des signes}}
\begin{prop}
Soit $P=a_nx^n+...+a_1x+a_0 \in \mathbb{R}[x]$, $r(P)$ le nombre de racines de
$P$ sur $\mathbb{R}^{*+}$, et $\sigma(P)$ le nombre de changements de signe
de la suite des coefficients $a_n,...,a_0$ de $P$. Alors $r(P) \leq \sigma(P)$.
\end{prop}
La preuve se fait par r\'ecurrence. Pour $n=0$ on a bien $0 \leq 0$.
Pour $n>0$ on applique l'hypoth\`ese de r\'ecurrence \`a $P'$ et on
utilise le fait que 1+le nombre de racines de $P'$ sur un intervalle
$]a,b[$ est sup\'erieur ou \'egal au nombre de racines de $P$ sur
$[a,b]$. Il y a deux cas possibles $\sigma(P')=\sigma(P)$ et
$\sigma(P')=\sigma(P)-1$. Dans le deuxi\`eme cas, comme $r(P')+1\geq
r(P)$, alors $r(P) \leq r(P')+1 \leq \sigma(P')+1=\sigma(P)$. Dans le
premier cas, il faut trouver une racine suppl\'ementaire pour
$P'$. Pour cela, on regarde ce qui se passe en $0^+$. Supposons que
$a_0>0$ pour fixer les id\'ees, comme $a_0$ est du m\^eme signe
que $a_1$ (ou que $a_2$ si $a_1=0$ etc.), $P'$ est positif en $0^+$
donc $P$ croit en $0^+$, donc doit atteindre un maximum local avant
la premi\`ere racine de $P$, ce maximum local est une racine de $P'$.

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

On d\'eduit de la r\`egle des signes de Descartes
un crit\`ere indiquant si un polyn\^ome poss\`ede 0, 1 ou 
plus de racines sur $]0,1[$ en se ramenant par changement de variables 
$x \leftarrow 1/x$ \`a $]1,+\infty[$ puis par $x \leftarrow x+1$
\`a $]0,+\infty[$. 

On prend alors le polyn\^ome de d\'epart et une majoration
sur les racines, par exemple $C=|P|_{\infty}/|\mbox{lcoeff}(P)|$. Si $z$
est une racine positive de $P$, alors $z/C$ est une racine positive
de $Q(x):=P(Cx)$, on est donc ramen\'e a chercher des intervalles
d'isolation de racines de $Q$ dans $]0,1[$. S'il y a 0 ou 1 changement
de signe, on conclut. Sinon, on teste si $1/2$ est racine
puis on cherche dans $]0,1/2[$ et dans $]1/2,1[$, ce qui
revient au probl\`eme pr\'ec\'edent sur $]0,1[$
avec les polyn\^omes $R(x)=Q(x/2)$ et $S(x)=Q((x+1)/2)$.

On fait de m\^eme sur $P(x=-x)$ pour les racines r\'eelles n\'egatives.
\begin{giaconload}
fonction descartes(P)
  local l,res,j;
  l:=coeffs(P);
  l:=remove(0,l);
  res:=0;
  pour j de 0 jusque dim(l)-2 faire
    si l[j]*l[j+1]<0 alors res++; fsi;
  fpour;
  return res;
ffonction:;
  
fonction isole01(P,a,b)
  // renvoie une liste d'intervalles d'isolation
  local n,m,Q,R,res;
  global x;
  Q:=x^degree(P)*P(x=1/x);
  R:=Q(x=x+1);
  n:=descartes(R);
  si n=0 alors return NULL; fsi;
  si n=1 alors return [a,b]; fsi;
  Q:=numer(P(x=x/2));
  R:=numer(P(x=(x+1)/2));
  m:=(a+b)/2;
  si P(x=1/2)=0 alors res:=(a+b)/2; sinon res:=NULL; fsi;
  res:=res,isole01(Q,a,m),isole01(R,m,b)
  retourne res;
ffonction:;

fonction isole(P) // racines positives de P
  local l,M;
  l:=coeffs(P);
  M:=ceil(maxnorm(l)/abs(l[0]))+1;
  P:=P(x=M*x);
  return isole01(P,0,M);
ffonction:;
\end{giaconload}
Racines de $P$ dans $\mathbb{R}^{*+}$
\giacinputmath{P:=(x-1/3)*(x-3)*(x^2+1); descartes(P)}
il y a un nombre pair de racines positives.
Racines de $P$ dans $]0,1[$
\giacinputmath{P1:=normal(x^degree(P)*P(x=1/x)); descartes(P1(x=x+1))}
Donc une racine dans $]0,1[$.

Isolation des racines positives de $P$: 
on commence par se ramener \`a chercher les racines
de $Q$ sur $]0,1[$ (il faudra les multiplier par $M$ \`a la fin)
\giacinputmath{M:=ceil(maxnorm(coeffs(P))/abs(lcoeff(P)))+1;Q:=P(x=M*x)}
\giacinputmath{P1:=normal(x^degree(Q)*Q(x=1/x)); descartes(P1(x=x+1))}
Il faut d\'ecouper en deux~:
\giacinputmath{Q1:=Q(x=x/2); P1:=normal(x^degree(Q1)*Q1(x=1/x)); descartes(P1(x=x+1))}
donc il y a une racine de $Q$ dans $]0,1/2[$, donc de $P$ dans $]0,5/2[$
\giacinputmath{Q1:=Q(x=(x+1)/2); P1:=normal(x^degree(Q1)*Q1(x=1/x)); descartes(P1(x=x+1))}
donc il y a une racine de $Q$ dans $]1/2,1[$, donc de $P$ dans $]5/2,5[$.
\giacinputmath{normal(Q(x=1/2))}
$Q$ n'est pas nul en $1/2$, on conclut donc qu'il existe 2 racines positives
pour $P$. V\'erification avec {\tt isole}~:
\giacinputmath{isole(P);}



\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 $\mathbb{R}$ et 4 de degr\'e 1
dans $\mathbb{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 \`a calculer les racines complexes approch\'ees
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 $\mathbb{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'\'egalit\'e est vraie. Soit 
\[ C=(f(x)-P(x))/\prod_{j=0}^n(x-x_j) \]
on consid\`ere 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 \`a $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\^eme 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\'e inf\'erieur ou \'egal \`a $n$ 
et $ \prod_{j=0}^n(x-x_j) - x^{n+1}$ est de degr\'e
inf\'erieur ou \'egal \`a $n$. Donc il existe bien un r\'eel $\xi_x$ dans
l'intervalle contenant les $x_j$ et $x$ tel que
\[ C=\frac{f^{[n+1]}(\xi_x)}{(n+1)!} \]

Illustration~: approcher la fonction $f(x)=x*\sin(x)-\cos(x)$ sur
$[-1,2]$ avec un pas de 0.5.\\
\giacinputmath{f(x):=x*sin(x)-cos(x);
X:=[-1.0,-0.5,0.0,0.5,1.0,1.5,2.0];
Y:=map(f,X);
P:=lagrange(X,Y);
}\\
\giacinput{plot([f(x),P],x=-2..3,color=[red,blue])}\\
Graphiquement, on ne distingue pas $f(x)$ et $P$ sur l'intervalle
d'interpolation. Tracons la diff\'erence~:
\giacinput{plot([f(x)-P],x=-1.2..2.2)}\\
D\`es qu'on sort de l'intervalle, l'erreur explose.
Comparons avec la majoration d'erreur sur l'intervalle\\
\giacinputmath{diff(f(x),x,7)}\\
La d\'eriv\'ee 7-i\`eme est donc major\'ee par 10 sur $[-1,2]$. On
trace la majoration de l'erreur avec l'erreur.
\giacinput{
Q:=10*product(x-X[j],j,0,size(X)-1);
plot([abs(f(x)-P),10*abs(Q)/7!],x=-1..2,color=[red,blue])
}\\
On observe qu'avec des points d'interpolation \'equir\'epartis,
l'erreur est plus grande aux bords de l'intervalle qu'au milieu.
Pour rendre la majoration de l'erreur uniforme, il faudrait utiliser
plus de points vers les bords de l'intervalle, c'est ce qu'on fait
sur $[-1,1]$ en utilisant les racines du polyn\^ome de Tchebycheff 
de 1\`ere esp\`ece $T_{n+1}$  d\'efini par
$$ T_k(\cos(x))=\cos(kx)$$
Ce polyn\^ome de degr\'e $k$ est uniform\'ement born\'ee par 1 sur $[-1,1]$.
Et si un polyn\^ome $U$ a le m\^eme degr\'e et coefficient dominant que $T_k$
alors sa norme  est au moins 1, car si elle \'etait strictement plus petite que 1,
la diff\'erence $T_k-U$ serait du signe de $T_k$ lorsque $T_k=\pm 1$
et changerait donc de signe $k$ fois sur l'intervalle $[-1,1]$ ce qui
est impossible pour un polyn\^ome de degr\'e au plus $k-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\'e pour exp.
%Programmation de la recherche de toutes les racines
%d'un polyn\^ome. Bassins d'attractions.

\section{Int\'egration num\'erique} % \'equations diff\'erentielles(?)
\index{integration} \index{quadrature}
Les fractions rationnelles admettent une primitive que l'on calcule
en d\'ecomposant la fraction avec B\'ezout comme expliqu\'e pr\'ec\'edemment.
Mais elles font figure d'exceptions,
la plupart des fonctions n'admettent pas de primitives qui s'expriment
\`a l'aide des fonctions usuelles. Pour calculer une int\'egrale,on
revient donc \`a la d\'efinition 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\'ecoupe l'intervalle d'int\'egration
en subdivisions $[a,b]=[a,a+h] + [a+h,a+2h]+...[a+(n-1)h,a+nh=b$, o\`u
$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\`ezes}
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 \`a droite et gauche)
ou $f((\alpha+\beta)/2)$ (point milieu), pour les trap\`ezes on utilise
le segment reliant $[\alpha,f(\alpha)]$ \`a $[\beta,f(\beta)]$.

Exemple~: calcul de la valeur approch\'ee de $\int_0^1 t^3 dt$
(on en connait la valeur exacte $1/4=0.25$) par ces m\'ethodes 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 \`a 9.
Pour les rectangles \`a 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 \`a 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\`ezes, l'aire du trap\`eze d\'elimit\'e par l'axe des $x$,
les verticales $y=\alpha$, $y=\beta$ et les points sur ces verticales
d'ordonn\'ees 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\`ezes, on voit que chaque terme apparait deux fois
sauf le premier et le dernier.

Plus g\'en\'erallement, 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\`u $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\`ezes donnent
une bien meilleure pr\'ecision que les rectangles. Plus g\'en\'eralement,
la pr\'ecision de l'approximation n'est pas la m\^eme selon le choix
de m\'ethode.
Ainsi pour les rectangles \`a gauche (le r\'esultat est le m\^eme
\`a droite), si $f$ est continument d\'erivable, de d\'eriv\'ee major\'ee
par une constante $M_1$ sur $[a,b]$, en faisant un
d\'eveloppement 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\'ee par $0.015$
sur une subdivision, donc par $0.15$ sur les 10 subdivisions.

Pour le point milieu, on fait le d\'eveloppement en $(\alpha+\beta)/2$ \`a l'ordre
2, en supposant que $f$ est deux fois continument d\'erivable~:
\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\'ee par $0.25e-3$, donc sur 10 subdivisions par $0.25e-2=0.0025$.

Pour les trap\`ezes, la fonction $g$ dont le graphe est le segment reliant
$[\alpha,f(\alpha)]$ \`a $[\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\'erivable, on peut donc majorer la diff\'erence
entre $f$ et $g$ en utilisant (\(\ref{eq:lagrange}\)), on int\`egre 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\`u $M_2$ est un majorant de $|f'{'}|$ sur $[a,b]$.

Lorsqu'on calcule l'int\'egrale sur $[a,b]$ par une de ces m\'ethodes,
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\'egrale~:
\begin{itemize}
\item pour les rectangles \`a 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\`ezes $M_2h^2 (b-a)/12$.
\end{itemize}
Lorsque $h$ tend vers 0, l'erreur tend vers 0, mais pas \`a la m\^eme vitesse,
plus rapidement pour les trap\`ezes et le point milieu
que pour les rectangles. Plus on
approche pr\'ecis\'ement la fonction sur une subdivision, plus la
puissance de $h$ va \^etre grande, plus la convergence sera rapide
lorsque $h$ sera petit, avec toutefois une contrainte fix\'ee par la
valeur de $M_k$, borne sur la d\'eriv\'ee $k$-i\`eme de $f$ (plus
$k$ est grand, plus $M_k$ est grand en g\'en\'eral). 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\'ethode} \index{ordre}
On appelle m\'ethode d'int\'egration l'\'ecriture d'une approximation
de l'int\'egrale 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\`u les $y_j$ sont dans l'intervalle $[\alpha,\beta]$, par exemple \'equir\'epartis
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\'ethode d'int\'egration est d'ordre $n$ si il y a \'egalit\'e
ci-dessus pour tous les polyn\^omes de degr\'e inf\'erieur ou \'egal \`a $n$
et non \'egalit\'e pour un polyn\^ome de degr\'e $n+1$.
Par exemple, les rectangles \`a 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\'ethode est d'ordre $n$ avec des $w_j\geq 0$ et 
si $f$ est $n+1$ fois continument d\'erivable,
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\'eveloppement de Taylor de $f$ par exemple
en $\alpha$ \`a 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\'ethode est exacte pour $T_n(f)$, on en d\'eduit 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\`es sommation sur les $n$ subdivisions, on obtient que~:
\begin{theorem}
Pour une m\'ethode d'ordre $n$ \`a 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\'eja trait\'es, mais pas forc\'ement 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\'ethode 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\`egre $t^3$ sur $[0,1]$ en 1 subdivision par cette m\'ethode,
on obtient
\[ \frac{1}{6} (0+ 4 \frac{1}{2^3} + 1)=\frac{1}{4} \]
c'est-\`a-dire le r\'esultat exact, ceci est aussi v\'erifi\'e pour $f$ polynome
de degr\'e inf\'erieur ou \'egal \`a 2 puisque l'approximation de Lagrange
de $f$ est alors \'egale \`a $f$. On en d\'eduit que la m\'ethode de Simpson
est d'ordre 3 (pas plus car la m\'ethode
de Simpson appliqu\'ee \`a l'int\'egrale de
$t^4$ sur $[0,1]$ n'est pas exacte). On peut m\^eme am\'eliorer (cf. par exemple
Demailly) la constante g\'en\'erale de la section pr\'ec\'edente pour la majoration
de l'erreur en~:
\[ |\int_a^b f - I(f)| \leq \frac{h^4}{2880} (b-a) M_4 \]
Cette m\'ethode n\'ecessite $2n+1$ \'evaluations de $f$ (le calcul
de $f$ est un point \'etant presque toujours 
l'op\'eration la plus couteuse en temps d'une
m\'ethode de quadrature), au lieu de $n$ pour les rectangles
et le point milieu et $n+1$ pour les trap\`ezes. Mais on a une majoration
en $h^4$ au lieu de $h^2$ donc le ``rapport qualit\'e-prix'' de la m\'ethode
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\'en\'eraliser l'id\'ee pr\'ec\'edente, d\'ecouper la subdivision 
$[\alpha,\beta]$ en $n$ parts \'egales et utiliser le polyn\^ome d'interpolation
en ces $n+1$ points $x_0=\alpha, x_1, ..., x_n=\beta$. 
Ce sont les m\'ethodes 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\`es 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\'ef\`ere utiliser la m\'ethode de Simpson en utilisant
un pas plus petit.

Il existe aussi d'autres m\'ethodes, par exemple les quadratures de Gauss
(on choisit d'interpoler en utilisant des points non \'equir\'epartis 
tels que l'ordre de la m\'ethode soit le plus grand possible)
ou la m\'ethode de Romberg qui est une m\'ethode d'acc\'el\'eration
de convergence bas\'ee sur la m\'ethode des trap\`ezes (on prend
la m\'ethode des trap\`ezes en 1 subdivision de $[a,b]$, puis
2, puis $2^2$, ..., et on \'elimine les puissances de $h$
du reste $\int f-I(f)$ en utilisant un th\'eor\`eme d'Euler-Mac Laurin
qui montre que le d\'eveloppement asymptotique de
l'erreur en fonction de $h$ ne contient que des puissances paires
de $h$). De plus, on peut \^etre amen\'e \`a faire varier le pas $h$
en fonction de la plus ou moins grande r\'egularit\'e 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:trapezes}\)) & 0 & 0 & $M_1 h (b-a)/2$ \\
point milieu & (\(\ref{eq:trapezes}\)) & 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\'erie/int\'egrale.
%S\'eries enti\`eres.
%Point milieu, RK.
%S\'eries enti\`eres. 


\section{Alg\`ebre lin\'eaire}
\subsection{Le pivot de Gauss}\index{pivot} \index{Gauss}
Cet algorithme permet de cr\'eer des z\'eros en effectuant des manipulations
r\'eversibles sur les lignes d'une matrice. Ces lignes peuvent
repr\'esenter les coefficients d'un syst\`eme lin\'eaire, on obtient
alors un syst\`eme lin\'eaire \'equivalent, ou les coordonn\'ees
d'un syst\`eme de vecteur, on obtient alors les coordonn\'ees
d'un syst\`eme de vecteur engendrant le m\^eme sous-espace vectoriel. On
peut \'egalement repr\'esenter 2 matrices $A$ et $B$ reli\'es par une relation
$Ax=B$, cette relation reste alors vraie au cours et donc
apr\`es la r\'eduction.

\subsubsection{L'algorithme}
L'algorithme est le suivant:
\begin{enumerate}
\item on initialise $c=1$ et $l=1$,
$c$ d\'esigne le num\'ero de colonne $c$ \`a r\'eduire, et $l$ le 
num\'ero de ligne \`a partir duquel on cherche un ``pivot''
(au d\'ebut $l$ et $c$ valent donc 1, en g\'en\'eral
les 2 augmentent de 1 \`a chaque it\'eration)
\item Si $c$ ou $l$ est plus grand que le nombre de colonnes ou
de lignes on arr\^ete.
\item  Si la colonne $c$ n'a
que des coefficients nuls \`a partir de la ligne $l$, on incr\'emente
le num\'ero de colonne $c$ de 1 et on passe \`a l'\'etape 2. Sinon,
on cherche la ligne dont le coefficient est en valeur absolue
le plus grand possible (en calcul approch\'e) ou le plus simple
possible (en calcul exact), on \'echange cette ligne avec la ligne $l$.
Puis on effectue pour toutes les lignes sauf $l$ ou pour toutes
les lignes \`a partir de $l+1$ (selon qu'il s'agit d'une r\'eduction
de Gauss compl\`ete ou d'une r\'eduction de Gauss sous-diagonale) la
manipulation r\'eversible 
\[ L_j \leftarrow  L_j -\frac{a_{jc}}{a_{lc}} L_l \]
On incr\'emente $c$ et $l$ de 1 et on passe \`a l'\'etape 2.
\end{enumerate}

\subsubsection{Efficacit\'e de l'algorithme}
Si la matrice poss\`ede $L$ lignes et $C$ colonnes,
le nombre maximal d'op\'erations pour r\'eduire une ligne est $C$ divisions,
$C$ multiplications, $C$ soustractions, donc $3C$ op\'erations
arithm\'etiques de base. Il y a $L-1$ lignes \`a r\'eduire \`a chaque
\'etape et min$(L,C)$ \'etapes \`a effectuer, on en d\'eduit que le
nombre maximal d'op\'erations pour r\'eduire une matrice est
$3LC$min$(L,C)$. Pour une matrice carr\'ee de taille $n$, cela fait 
$3n^3$ op\'erations.

\subsubsection{Erreurs d'arrondis du pivot de Gauss} \index{erreur}
Comme $|a_{jc}| \leq |a_{lc}|$, une \'etape de r\'eduction multiplie
au plus l'erreur absolue des coefficients par 2. Donc la
r\'eduction compl\`ete d'une matrice peut multiplier au pire l'erreur
absolue sur les coefficients par $2^n$ (o\`u $n$ est le nombre
d'\'etapes de r\'eduction, inf\'erieur au plus petit du nombre de lignes
et de colonnes). Ceci signifie qu'avec la pr\'ecision d'un double,
on peut au pire perdre toute pr\'ecision pour des matrices pas
si grandes que \c{c}a ($n=52$). Heureusement, il semble qu'en pratique, 
l'erreur absolue
ne soit que tr\`es rarement multipli\'ee par un facteur sup\'erieur \`a 10.

Par contre, si on ne prend pas la pr\'ecaution de choisir le pivot
de norme maximale dans la colonne, les erreurs d'arrondis se
comportent de mani\`ere 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\'eresser \`a la sensibilit\'e de la solution d'un
syst\`eme lin\'eaire \`a des variations de son second membre. Le traitement
du sujet \`a ce niveau est un peu difficile, cela fait intervenir
le nombre de conditionnement de la matrice $A$ du syst\`eme (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\'ecision).

\subsection{Applications de Gauss}
\subsubsection{Base d'un sous-espace} 
On r\'eduit la matrice des vecteurs \'ecrits en ligne, puis
on prend les lignes non nulles, dont les vecteurs
forment une base du sous-espace vectoriel engendr\'e par
les vecteurs du d\'epart.

Exemple~:  base du sous-espac engendr\'e par $(1,2,3), (4,5,6), (7,8,9)$.
On r\'eduit la matrice, la 3\`eme ligne est nulle
donc on ne garde que les 2 premi\`eres lignes $(1,0,-1), (0,1,2)$
(remarque: une r\'eduction sous-diagonale aurait suffi).

\subsubsection{D\'eterminant} \index{determinant}
On r\'eduit la matrice (carr\'ee) 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\'eduction sous forme \'echelonn\'ee (rref)} \index{rref} 
\index{reduction}
On r\'eduit la matrice puis on divise chaque ligne par son premier
coefficient non nul. Si la matrice repr\'esentait un syst\`eme
lin\'eaire inversible on obtient la matrice identit\'e sur les colonnes
sauf la derni\`ere et la solution en lisant la derni\`ere colonne.
La relation conserv\'ee est en effet $Ax=b$ o\`u $x$ est la solution
de l'\'equation, et \`a la fin de la r\'eduction $A=I$.

Par exemple pour r\'esoudre le syst\`eme
\[ \left\{ \begin{matrix}{ccc}
  x + 2y + 3z &=& 5 \\
  4x + 5y + 6z &=& 0 \\
  7x + 8y &=& 1
\end{matrix}
\right. \]
on r\'eduit sous forme \'echelonn\'ee 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\`u la solution
$x=-9, y=8, z=-2/3$.

\subsubsection{Inverse} \index{inverse}
On accolle la matrice identit\'e \`a droite de la matrice \`a inverser. On effectue
la r\'eduction sous forme \'echelonn\'ee, on doit obtenir \`a droite
l'identit\'e si la matrice est inversible, on a alors \`a gauche
la matrice inverse.
La relation conserv\'ee est en effet $A x=B$  o\`u $x$ est l'inverse
de la matrice de d\'epart, et en fin de r\'eduction $A=I$.

Par exemple, pour calculer l'inverse
de \verb|[[1,2,3],[4,5,6],[7,8,0]]|, on r\'eduit 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\'eduit la matrice sous forme \'echelonn\'ee. Puis on introduit des lignes
de 0 afin que les 1 en t\^ete de ligne se trouvent sur la diagonale
de la matrice. On enl\`eve ou on rajoute des lignes de 0 \`a la fin
pour obtenir une matrice carr\'ee. Une base du noyau est alors
form\'ee en prenant chaque colonne correspondant \`a un 0 sur la
diagonale, en rempla\c{c}ant ce 0 par -1. On v\'erifie qu'on obtient bien 0
en faisant le produit de ce vecteur par la matrice r\'eduite. De
plus les vecteurs cr\'e\'es sont clairement lin\'eairement ind\'ependants
(puisqu'ils sont \'echelonn\'es), et il y en a le bon nombre (th\'eor\`eme
noyau-image).

Exemple~: calcul du noyau de \verb|[[1,2,3,4],[1,2,7,12]]|, on
r\'eduit 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\`eme 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 \`a la fin pour rendre la matrice carr\'ee.
On obtient ainsi le syst\`eme \'equivalent de matrice
\verb|[[1,2,0,-2],[0,0,0,0],[0,0,1,2],[0,0,0,0]]|.
La 2\`eme colonne donne le premier vecteur de la base du noyau,
$(2,-1,0,0)$, la 4\`eme colonne donne le deuxi\`eme vecteur
$(-2,0,2,-1)$, on v\'erifie ais\'ement que ces 2 vecteurs
forment une famille libre du noyau, donc une base car 
la dimension du noyau est \'egale \`a 2 (dimension de l'espace de d\'epart
moins le rang de la matrice, c'est-\`a-dire le nombre de lignes 
non nulles de la matrice r\'eduite).

\subsubsection{La m\'ethode de factorisation $LU$} \index{LU}
Nous ne la d\'eveloppons pas \`a ce niveau, elle permet d'\'ecrire
une matrice $A$ comme produit de deux matrices triangulaire
inf\'erieures et sup\'erieures, ce qui permet de ramener la r\'esolution
de syst\`eme \`a la r\'esolution de deux syst\`emes triangulaires.

\subsection{R\'eduction exacte des endomorphismes} \index{diagonalisation}
On calcule le polynome caract\'eristique 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\'ethodes \'evitant
le calcul de noyau, m\'ethode de Fadeev-Laguerre-Souriau que
nous ne pr\'esentons pas ici.

\subsubsection{Polynome caract\'eristique} \index{polynome caracteristique}
On peut le calculer en d\'eveloppant le d\'eterminant det$A-\lambda I$, mais
il est plus efficace de le calculer par interpolation.
Soit $A$ une matrice carr\'ee de taille $n$, 
on sait que son polynome caract\'eristique
est un polynome de degr\'e $n$, il suffit
de connaitre sa valeur en $n+1$ points distincts, on calcule
donc $n+1$ d\'eterminants det$A-\lambda I$ en rempla\c{c}ant $\lambda$ par
sa valeur (il y a plus de d\'eterminants \`a calculer mais
ce sont des d\'eterminants sans param\`etre $\lambda$ donc beaucoup
plus simple \`a calculer), ce qui permet de reconstruire le polynome
caract\'eristique 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\'erifier que c'est bien le polynome caract\'eristique.

Il faut effectuer $n+1$ calculs de d\'eterminants, ce qui n\'ecessite
$O(n^4)$ op\'erations. Il existe des m\'ethodes plus efficaces, par
exemple le calcul du polynome minimal probabiliste pr\'esent\'e
plus bas ($O(n^3)$ op\'erations).

\subsubsection{Polynome minimal} \index{polynome minimal}
\begin{defn}
Le polynome minimal d'une matrice $A$ est un polyn\^ome $M$ de degr\'e minimal
tel que $M(A)=0$ et de coefficient dominant \'egal \`a 1. 
Un tel polynome divise tous les polynomes
tels que $P(A)=0$, il divise le polynome caract\'eristique
de $A$ et il a les m\^emes racines que le polynome caract\'eristique.
\end{defn}
Preuve:\\
D'abord $M$ divise tous les polynomes
tels que $P(A)=0$, car si $R$ d\'esigne 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\'e est plus petit que celui de $M$.

En particulier le polynome minimal divise le polynome caract\'eristique $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\'eterminant fois l'identit\'e, soit $C(\lambda)I$. Comme
$C(\lambda)I-C(A)$ peut se factoriser par $\lambda I -A$ en appliquant
(\(\ref{eq:factor}\)) \`a chaque monome de $C$, on
en d\'eduit que $C(A)$ se factorise par $\lambda I -A$, donc $C(A)=0$
en regardant les termes de plus haut degr\'e de ces polynomes en
$\lambda$ \`a coefficients matriciels).

Montrons enfin que les racines du polynome caract\'eristique sont racines
du polynome minimal. En effet soit $\lambda$ une racine
du polynome caract\'eristique 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 \^etre inversible. Comme $M(A)=0$
on en d\'eduit que $M(\lambda)I$ n'est pas inversible donc $M(\lambda)=0$,
$\lambda$ est une racine de $M$.
Donc si le polynome caract\'eristique n'a pas de racines multiples,
il est \'egal au polynome minimal.

Pour calculer $M$, on peut chercher une relation de degr\'e minimal entre
les puissances de $A$, en les voyant comme des vecteurs \`a $n^2$
composantes (ce qui revient \`a aplatir en un long vecteur
tous les coefficients de la matrice). Cela revient \`a calculer
le noyau de l'application lin\'eaire dont les colonnes sont les coefficients
des puissances de $A$ (de 0 \`a $n$), en gardant le premier vecteur
obtenu par l'algorithme calcul du noyau ci-dessus.

Cette m\'ethode est toutefois couteuse, car il faut r\'eduire une matrice
ayant $n^2$ lignes et $n+1$ colonnes. Il existe une autre m\'ethode
moins couteuse et qui marche presque toujours. Elle consiste \`a
calculer le polynome minimal de $A$ par rapport \`a un vecteur $v$
c'est-\`a-dire le polynome de degr\'e 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\'eristique. Si par
chance, on trouve que $M_v$ est de degr\'e $n$, alors $M_v$ sera
\'egal \`a $M$ et au polynome caract\'eristique. On fait donc le calcul
du noyau de l'application lin\'eaire dont les colonnes sont les $A^jv$
pour $j$ variant de 0 \`a $n$. Si l'on trouve un espace de dimension 1,
alors $M_v$ est de degr\'e $n$ et on a simultan\'ement le polynome
minimal et caract\'eristique avec le polynome correspondant \`a 
ce vecteur du noyau. Si le degr\'e 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\'e $n$
on conclut, sinon on peut tester si ce polynome \'evalu\'e en $A$
est nul, ce sera alors le polynome minimal.

Exemple, on reprend la matrice \verb|[[1,-1],[2,4]]|,
et comme vecteur al\'eatoire $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 \'ecrit 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\`es) $-6+5x-x^2$. 
Comme il est de degr\'e maximal 2, c'est le polynome minimal
et caract\'eristique.

\subsection{R\'eduction approch\'ee des endomorphismes}
On pourrait trouver des valeurs propres approch\'ees d'une matrice
en calculant le polynome caract\'eristique ou minimal puis en le
factorisant num\'eriquement. Mais cette m\'ethode n'est pas id\'eale
relativement aux erreurs d'arrondis (calcul du polynome caract\'eristiaue,
de ses racines, et nouvelle approximation en calculant le noyau
de $A-\lambda I$), lorsqu'on veut calculer quelques valeurs propres
on pr\'ef\`ere utiliser des m\'ethodes it\'eratives directement sur $A$
ce qui \'evite la propagation des erreurs d'arrondi.

\subsubsection{M\'ethode de la puissance} \index{puissance}
Elle permet de d\'eterminer 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\'eatoire $v$ et on calcule la suite
$v_n=Av_{n-1}=A^n v $. Si $v$ a pour coordonn\'ees $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\`ese 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\'eom\'etrique de raison $x_j/x_n$ converge vers 0.
Autrement dit, si $V_n\neq 0$ (ce qui a une probabilit\'e 1 d'\^etre
vrai pour un vecteur al\'eatoire), 
$v_n$ est \'equivalent \`a $V_n x_n^n e_n$. Lorsque $n$
est grand, $v_n$ est presque colin\'eaire au vecteur propre
$e_n$ (que l'on peut prendre comme $v_n$ divis\'e par sa norme), 
ce que l'on d\'etecte en testant si $v_{n+1}$ et $v_n$
sont presques colin\'eaires, et de plus le facteur de colin\'earit\'e
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\'e \`a $\lambda=3$ pour la matrice \verb|[[1,-1],[2,4]]|
et le vecteur $(1,0)$.

Lorsqu'on applique cette m\'ethode a une matrice r\'eelle, il peut arriver
qu'il y ait deux valeurs propres conjugu\'ees de module maximal. Le m\^eme
type de raisonnement montre que pour $n$ grand,
$v_{n+2}$ est presque colin\'eaire \`a l'espace engendr\'e 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\'erie g\'eom\'etrique, on gagne le m\^eme nombre
de d\'ecimales \`a chaque it\'eration.

\subsubsection{It\'erations inverses} \index{iterations inverses}
La m\'ethode pr\'ec\'edente permet de calculer la valeur propre
de module maximal d'une matrice. Pour trouver une valeur
propre proche d'une quantit\'e donn\'ee $x$, on peut appliquer
la m\'ethode pr\'ec\'edente \`a 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\'ees}
Si la matrice $A$ est sym\'etrique, et
si $e_n$ est un vecteur propre norm\'e \'ecrit en colonne, on peut consid\'erer
la matrice $B=A-x_n e_n e_n^t$ qui poss\`ede les m\^emes valeurs
propres et m\^emes vecteurs propres que $A$ avec m\^eme multiplicit\'e, 
sauf $x_n$ qui est remplac\'e 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\`eme valeur propre (en valeur absolue), l'\'eliminer
et ainsi de suite.

Si la matrice $A$ n'est pas sym\'etrique, on peut utiliser une technique
analogue lorsque 0 n'est pas valeur propre de $A$ (on peut s'y
ramener en ajoutant \`a $A$ un multiple de l'identit\'e).
En effet on peut construire un vecteur propre de $B$ pour une
valeur propre $x_k \neq 0$ \`a 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\`u l'\'equation
\[ y x_k  = x_n e_n.e_k \]
N\'eanmoins cette m\'ethode 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{Guide rapide KhiCAS sur calculatrices}
\subsection{Installation}
\subsubsection{Installation Casio Graph 90}
(Ceci s'applique dans une certaine mesure aux Casio Graph 35eii).
Rendez-vous ici
\ahref{https://www-fourier.univ-grenoble-alpes.fr/\home{parisse}/install_fr.html#ti}{pour
installer}.
\ahref{https://www-fourier.univ-grenoble-alpes.fr/\home{parisse}/casio/khicasio.html}{Documentation
  Casio}
Si Xcas crashe sur la Casio, essayez d'abord
de taper sur la touche MENU puis lancez une autre
application puis à nouveau MENU puis relancez Xcas. Si cela
ne fonctionne pas ou si Xcas est bloqué, redémarrer la calculatrice
(éventuellement en appuyant avec une pointe sur le bouton RESET au dos).
Si Xcas bloque au lancement, allez dans le MENU principal, sélectionnez
l'application Memoire, puis mémoire de stockage, sélectionnez
le fichier \verb|session.xw| et effacez-le.

\subsubsection{Installation TI Nspire CX/CX II}
Vous pouvez installer KhiCAS sur TI Nspire CX/CX II \`a condition de
ne pas avoir install\'e la toute derni\`ere mise \`a jour de l'OS de TI. Rendez-vous ici
\ahref{https://www-fourier.univ-grenoble-alpes.fr/\home{parisse}/install_fr.html#ti}{pour
installer}.
\ahref{https://www-fourier.univ-grenoble-alpes.fr/\home{parisse}/ti/khicasti.html}{Documentation
Nspire}.

\subsubsection{Installation Numworks N0110 non verrouill\'ee.}
Pour les Numworks, il faut un mod\`ele N0110 (juillet 2019 ou plus
tard) et ne surtout pas avoir fait de mise \`a jour vers Epsilon 16 ou
sup\'erieur. Rendez-vous sur
\ahref{https://www-fourier.univ-grenoble-alpes.fr/\home{parisse}/nws.html}{cette
  page} avec un navigateur compatible web-USB.
\ahref{https://www-fourier.univ-grenoble-alpes.fr/\home{parisse}/numworks/khicasnw.html}{Documentation Numworks}.

\subsection{Le shell KhiCAS}
Si KhiCAS n'est pas lanc\'e~:
\begin{itemize}
\item Sur Casio, depuis le menu (touche MENU), s\'electionner
Xcas puis EXE.
\item Sur TI Nspire,
touche maison, puis Parcourir, puis depuis le r\'epertoire ndless,
s\'electionner khicas.tns et valider. Si une erreur de format non
pris en charge apparait, (re-)lancez l'installation de ndless.
\item Sur Numworks, touche Maison puis OK.
\end{itemize}
Pour \'eteindre la calculatrice, il faut quitter KhiCAS sur les TI
Nspire (taper doc deux fois avant ctrl on) et les Casio (taper MENU
avant shift AC/ON).

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

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

Pour bloquer le clavier en mode alphab\'etique~: sur les Casio, taper
sur la touche F5, sur les Numworks taper 2 fois la touche alpha.
Pour quitter le mode alphab\'etique, taper sur la touche alpha.

Exemples~: 
\begin{enumerate}
\item pour d\'efinir une fonction, sur Casio faire shift F6 (ou
shift PRGM), sur TI/Numworks faire shift ),
puis s\'electionner \verb|f(x):=| et entrer l'expression
de la fonction. Changez si n\'ecessaire le nom de fonction.
\item pour d\'efinir une expression (attention \`a la diff\'erence
  avec une fonction), taper l'expression puis
sur Casio la touche fl\`eche vers
la droite (au-dessus de AC/ON) puis F5 puis le nom de variable,
sur TI ctrl sto, sur Numworks shift sto.
\item Pour d\'eriver une fonction ou une expression par rapport \`a
  $x$,  on peut utiliser \verb|'| (menu F2
  ou shift-2). Sinon utiliser \verb|diff| (m\^eme menu),
  par exemple \verb|diff(f,x,4)| pour avoir la dérivée 4-ième
  d'une expression.
\item  Par exemple pour obtenir les 4 premiers termes de
  la méthode de Newton pour résoudre $f(x)=0$ avec $u_0=2.3$,
  entrer la fonction $f$ (cf. 1 ci-dessus), puis $g=f'$ (cf 3 ci-dessus),
  puis taper \verb|U:=2.3|, puis taper \verb|U:=U-f(U)/g(U)|,
  puis recopier dans l'historique cette commande et exécutez la
  3 fois.
\item Pour tracer le graphe d'une fonction, saisir la commande
  \verb|plot(| depuis shift-F4 (casio) ou shift-3 (TI/Numworks).
Tapez l'expression \`a tracer. 
Ou fl\`eche vers le bas pour voir de l'aide et des exemples.
\item Les commandes les plus utiles sont accessibles par les
menus rapides via les raccourcis claviers. 
Pour voir plus de commandes~: F4 sur Casio, menu sur Nspire, touche
Toolbox sur Numworks.
\item Pour entrer ou \'editer une matrice,
  taper shift Matr sur Casio ou shift 7 sur
  TI/Numworks puis 1 matrice, donner le nom de la matrice, puis
nombre de lignes/colonnes pour une nouvelle matrice, puis entrer
les coefficients un par un puis EXE ou OK ou enter.
\end{enumerate}

\subsection{Programmation}
Pour saisir un programme, depuis le shell KhiCAS, taper sur EXIT (Casio) ou esc (Nspire) ou
Back (Numworks). Vous pouvez choisir entre syntaxe Xcas ou Python
(shift SETUP sur les Casio). Le parser d\'etecte aussi automatiquement
la syntaxe, par exemple si vous commencez \`a taper \verb|def | en
d\'ebut de programme. Les menus rapides F1 \`a F3 sur Casio et shift 1 \`a
shift 3 sur TI/Numworks contiennent les structures de programmation
les plus utiles.

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

Puis vous pouvez taper EXIT, esc ou Back pour
revenir au shell, d'o\`u vous pouvez ex\'ecuter le programme. 
On peut interrompre l'ex\'ecution d'un programme en appuyant sur
la touche AC/ON ou ON ou Back selon la calculatrice.
Si un programme ne fonctionne pas comme pr\'evu, la commande
\verb|debug()| permet de l'ex\'ecuter en pas \`a pas en affichant
les variables locales et permet de le mettre au point.

\subsection{Compatibilité avec Xcas}
Casio Graph 90 et 35:
Vous pouvez sauvegarder depuis Xcas une session de travail par le menu
Fich puis Exporter comme, puis Khicas Casio. On
peut ensuite transférer le fichier vers la calculatrice, en la
branchant au PC avec un cable prévu pour calculatrices Casio graph ou TI83,
(la calculatrice apparait alors comme un périphérique de stockage USB).
Puis dans KhiCAS sur la calculatrice, ouvrir la session. Dans l'autre
sens, transférer la session depuis la calculatrice et ouvrez-là sur le PC.
Lors de la sauvegarde sur calculatrice, pour des sessions pas trop longues,
un QR code est affiché, permettant de continuer les calculs sur
la version web de Xcas.

TI Nspire: Vous pouvez sauvegarder depuis Xcas une session de travail
par le menu Fich puis Exporter comme, puis session Khicas TI. On
peut ensuite transférer le fichier vers la calculatrice avec le logiciel
de connectivité TI ou avec le logiciel libre n-link. Dans l'autre
sens, transférer la session depuis la calculatrice et ouvrez-là sur le PC.

Numworks: on peut directement transférer une session depuis le menu
Fich, Numworks si la calculatrice est connectée.

\section{Quelques r\'ef\'erences}
\begin{itemize}
\item Analyse num\'erique et \'equations diff\'erentielles, 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\'ematiques concr\`etes, illustr\'ees 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\'en\'eralit\'e puisque \'echanger $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\`es (\(\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\'ecroissante 
minor\'ee (par $v_0$), $(v_n)$ est croissante major\'ee (par $u_0$), ces 2 suites sont 
convergentes et comme $u_{n+1}=\frac{u_n+v_n}{2}$, elles convergent vers la m\^eme 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\'ecisons maintenant la vitesse de convergence lorsque $a \geq b>0$. 
On va commencer par estimer le nombre
d'it\'erations n\'ecessaires pour que $u_n$ et $v_n$ soient du m\^eme 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\'erations $n_0$ est proportionnel
au log du log du rapport $a/b$.
Ensuite on est ramen\'e \`a \'etudier la convergence de la suite arithm\'etico-g\'eom\'etrique
de premiers termes $a=u_{n_0}$ et $b=v_{n_0}$ et m\^eme en tenant compte
de $M(a,b)=aM(1,b/a)$ \`a $a=1$ et $b=v_n/u_n$ donc $0\leq a-b \leq 1-e^{-0.1}$.
Alors l'\'equation (\(\ref{eq:agm1}\)) entraine 
\[ u_{n+1}-v_{n+1} \leq \frac{1}{8}(u_n-v_n)^2 \]
puis (par r\'ecurrence)
\[ 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\'erieure \`a une pr\'ecision donn\'ee $\epsilon$
au bout d'un nombre d'it\'erations 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\'eterminant $n$ chiffres significatifs,
il faudra alors $O(\ln(n))$ it\'erations pour se ramener \`a $M(1,b)$ avec $b\in [e^{-0.1},1]$ 
puis $O(\ln(n))$ it\'erations pour avoir la limite avec $n$ chiffres significatifs.

{\bf Le cas complexe}\\
On suppose maintenant que $a, b \in \mathbb{C}$ avec $\Re(a)>0, \Re(b)>0$. On va voir que
la suite arithm\'etico-g\'eom\'etrique converge encore. \\
{\bf \'Etude de l'argument}\\
On voit ais\'ement (par r\'ecurrence)
que $\Re(u_n)>0$ ; de plus $\Re(v_n) > 0$ car par d\'efinition de la racine carr\'ee
$\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 \^etre un r\'eel n\'egatif.
On en d\'eduit 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\`es $n$ it\'erations, on a 
\[ |\arg(u_n)-\arg(v_n)| \leq \frac{\pi}{2^n} \]
Apr\`es quelques it\'erations, $u_n$ et $v_n$ seront donc presque align\'es. 
Faisons 4 it\'erations.
On peut factoriser par exemple $v_n$ et on
est ramen\'e \`a l'\'etude de la suite de termes initiaux $a=u_n/v_n$ d'argument 
$\arg(u_n)-\arg(v_n)$ petit 
(inf\'erieur en valeur absolue \`a $\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\'esigne 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\'eduit
\[ |\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\'eom\'etrique, 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\'efinie si $a,b \in \mathbb{C}$ v\'erifient $\Re(a)>0, \Re(b)>0$,
on peut montrer que la relation ci-dessus s'\'etend (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\`u 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\'ecompose l'int\'egrale en 2 parties $[0,k]$ et $[k,\pi/2]$.
Sur $[0,k]$ on utilise (\(\ref{eq:maj1}\)), on majore chaque terme s\'epar\'ement
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\`u~:
\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\'egration 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\`eme 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\`u $J$ est donn\'e 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\'erieur \`a 1 et sup\'erieur \`a
\[ 
\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\'erieur \`a
\[ 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\'e\'ecrire
\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\'erations)
$M(1,2^-{27})=M_1=0.0781441403763$. 
On a, avec une erreur inf\'erieure \`a $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\'ee de $\ln(2)$ et r\'eciproquement.
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{giacjshere}
\end{document}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: 
