% Jean-Pierre Demailly
% Universit\'e de Grenoble I, Institut Fourier
% \'Equations Diff\'erentielles et Analyse Num\'erique
% Collection Grenoble-Sciences, Presses Universitaires de Grenoble
% Chapitre I

\input equa-mac.tex

\blankline\blankline
\hfill\chtitle Chapitre I|%
\blankline\blankline
\hfill\title Calculs num\'eriques approch\'es|
\chaptitlerunning I -- Calculs num\'eriques approch\'es|
\vskip4pc

L'objet de ce chapitre est de mettre en \'evidence les principales
difficult\'es li\'ees \`a la pratique des calculs num\'eriques sur
ordinateur. Dans beaucoup de situations, il existe des m\'ethodes
sp\'ecifiques permettant d'accro\^{\i}tre \`a la fois l'efficacit\'e et
la pr\'ecision des calculs.

\section 1. Cumulation des erreurs d'arrondi|

\subsection 1.1. Repr\'esentation d\'ecimale approch\'ee des
nombres r\'eels|

\noindent La capacit\'e m\'emoire d'un ordinateur est par construction
finie. Il est donc n\'ecessaire de repr\'esenter les nombres r\'eels 
sous forme approch\'ee. La notation la plus utilis\'ee \`a l'heure actuelle
est la repr\'esentation avec virgule flottante~: un nombre r\'eel $x$ est
cod\'e sous la forme
$$
x\simeq \pm\kern2pt m\cdot b^p
$$
o\`u $b$ est la {\it base de num\'eration}, $m$ la {\it mantisse}, et
$p$ l'exposant. Les calculs internes sont g\'en\'eralement effectu\'es
en base $b=2$, m\^eme si les r\'esultats affich\'es sont finalement
traduits en base 10.

La mantisse $m$ est un  nombre \'ecrit avec virgule fixe et poss\'edant
un nombre maximum $N$ de chiffres significatifs (impos\'e par le choix
de la taille des emplacements m\'emoires allou\'es au type {\it
r\'eel})~: suivant les machines, $m$ s'\'ecrira
\vskip9pt
$\matrix{
\displaystyle \bullet\quad m=0,a_1a_2\dots a_N=\sum_{k=1}^N a_kb^{-k},\hfill
&b^{-1}\leq m<1~;\hfill\cr
\displaystyle \bullet\quad m=a_0,a_1a_2\dots a_{N-1}=\sum_{0\leq k<N}
a_kb^{-k},\hfill&1\leq m <b.\hfill\cr
}$
\vskip9pt
\noindent Ceci entra\^{\i}ne que la pr\'ecision dans l'approximation
d'un nombre r\'eel est toujours une {\it pr\'ecision relative}~:

\bigskip
\majdeb
$$
\dfrac{\Delta x}{x}=\dfrac{\Delta m}{m}\leq\dfrac{b^{-N}}{b^{-1}}=b^{1-N}.
$$
\vskip-\baselineskip
\majfin
\vskip\baselineskip

\noindent On notera $\varepsilon=b^{1-N}$ cette pr\'ecision relative.

En Turbo Pascal, les r\'eels occupent 6 octets de m\'emoire. Ceci
permet de repr\'esenter les r\'eels avec une mantisse de 11 chiffres
significatifs, dans une \'echelle allant de $10^{-38}=1\ {\rm E}-38$
\`a $10^{38}=1\ {\rm E}+38$. La pr\'ecision relative est donc de
l'ordre de $10^{-10}$.

\subsection 1.2. Non-associativit\'e des op\'erations arithm\'etiques|

\noindent Supposons par exemple que les r\'eels soient calcul\'ees avec
3 chiffres significatifs et arrondis \`a la d\'ecimale la plus proche.
Soit \`a calculer la somme $x+y+z$ avec
$$
x=8,22,\quad y=0,00317,\quad z=0,00432
$$
$\matrix{
(x+y)+z\ {\rm donne}~:&x+y=8,22317\simeq 8,22\hfill\cr
\noalign{\vskip3pt}
&(x+y)+z\simeq 8,22432\simeq 8,22\hfill\cr
\noalign{\vskip3pt}
x+(y+z)\ {\rm donne}~:&y+z=0,00749\hfill\cr
\noalign{\vskip3pt}
&x+(y+z)=8,22749\simeq 8,23.\hfill\cr
}$
\vskip6pt
\noindent L'addition est donc non associative par suite des erreurs d'arrondi~!

\subsection 1.3. Erreur d'arrondi sur une somme|

\noindent On se propose d'\'etudier quelques m\'ethodes permettant de {\it
majorer} les
erreurs d'arrondi dues aux op\'erations arithm\'etiques.

Soient $x,y$ des nombres r\'eels suppos\'es repr\'esent\'es sans erreur avec
$N$
chiffres significatifs~:
$$
\matrix{
x=0,a_1a_2\dots a_N\cdot b^p,\hfill&b^{-1+p}\leq x<b^p\hfill\cr
\noalign{\vskip3pt}
y=0,a'_1a'_2\dots a'_N\cdot b^q,\hfill&b^{-1+q}\leq y<b^q\hfill\cr
}
$$
Notons $\Delta(x+y)$ l'erreur d'arrondi commise sur le calcul de $x+y$.
Supposons par
exemple $p\geq q$. S'il n'y a pas d\'ebordement, c'est-\`a-dire si $x+y<b^p$,
le calcul
de $x+y$ s'accompagne d'une perte des $p-q$ derniers chiffres de $y$
correspondant aux
puissances $b^{-k+q}<b^{-N+p}$~; donc $\Delta(x+y)\leq b^{-N+p}$, alors que
$x+y\geq
x\geq b^{-1+p}$. En cas de d\'ebordement $x+y\geq b^p$ (ce qui se produit par
exemple si
$p=q$ et $a_1+a'_1\geq b$), la d\'ecimale correspondant \`a la puissance
$b^{-N+p}$ est
elle aussi perdue, d'o\`u $\Delta(x+y)\leq b^{1-N+p}$. Dans les deux cas~:

\bigskip
\majdeb
\centerline{$
\Delta(x+y)\leq\varepsilon(|x|+|y|)
$,}
\par\majfin\par

\bigskip
\noindent  o\`u $\varepsilon=b^{1-N}$ est la
pr\'ecision relative d\'ecrite au \S 1.1. Ceci reste vrai quel que soit le
signe des nombres
$x$ et $y$.

En g\'en\'eral, les r\'eels $x,y$ ne sont  eux-m\^emes connus que par des
valeurs
approch\'ees $x',y'$ avec des erreurs respectives $\Delta x=|x'-x|$, $\Delta
y=|y'-y|$. A
ces erreurs s'ajoute l'erreur d'arrondi
$$
\Delta(x'+y')\leq \varepsilon(|x'|+|y'|)\leq \varepsilon(|x|+|y|+\Delta
x+\Delta y).
$$
Les erreurs $\Delta x$, $\Delta y$ sont elles-m\^emes le plus souvent d'ordre
$\varepsilon$ par rapport \`a $|x|$ et $|y|$, de sorte que l'on pourra
n\'egliger les
termes $\varepsilon\Delta x$ et $\varepsilon\Delta y$. On aura donc~:

\bigskip
\majdeb
\centerline{$
\Delta(x+y)\leq \Delta x+\Delta y+\varepsilon(|x|+|y|).
$}
\par\majfin\par

\bigskip
\noindent Soit plus g\'en\'eralement \`a calculer une somme
$\displaystyle\sum_{k=1}^nu_k$
de r\'eels {\it positifs}. Les sommes partielles $s_k=u_1+u_2+\dots+u_k$ vont
se
calculer de proche en proche par les formules de r\'ecurrence
$$
\cases{
s_0=0\hfill\cr
\noalign{\vskip3pt}
s_k=s_{k-1}+u_k,\hfill&\quad$k\geq 1$.\hfill\cr
}
$$
Si les r\'eels $u_k$ sont connus exactement, on aura sur les sommes $s_k$ des
erreurs
$\Delta s_k$ telles que $\Delta s_1=0$ et
$$
\Delta s_k\leq \Delta s_{k-1}+\varepsilon(s_{k-1}+u_k)=\Delta
s_{k-1}+\varepsilon s_k.
$$
L'erreur globale sur $s_n$ v\'erifie donc
$$
\Delta s_n\leq \varepsilon(s_2+s_3+\dots+s_n),
$$
soit
$$
\Delta s_n\leq \varepsilon(u_n+2u_{n-1}+3u_{n-2}+\dots+(n-1)u_2+(n-1)u_1).
$$
Comme ce sont les premiers termes somm\'es qui sont affect\'es des plus gros
coefficients dans l'erreur $\Delta s_n$, on en d\'eduit la r\`egle g\'en\'erale
suivante
(cf. exemple 1.2).

\bigskip
\majdeb
\state R\`egle g\'en\'erale -- | Dans une sommation de r\'eels,
l'erreur a tendance \`a \^etre minimis\'ee lorsqu'on somme en premier
les termes ayant la plus petite valeur absolue.
\par\majfin
\par

\subsection 1.4. Erreur d'arrondi sur un produit|

\noindent Le produit de deux mantisses de $N$ chiffres donne une mantisse de
$2N$ ou $2N-1$ chiffres dont les $N$ ou $N-1$ derniers vont \^etre
perdus. Dans le calcul d'un produit $xy$ (o\`u $x,y$ sont suppos\'es
repr\'esent\'es sans erreur) il y aura donc une erreur d'arrondi
$$
\Delta(xy)\leq \varepsilon|xy|,\quad\hbox{o\`u}\quad \varepsilon=b^{1-N}.
$$
Si $x$ et $y$ ne sont eux-m\^emes connus que par des valeurs approch\'ees $x'$,
$y'$ et
si $\Delta x=|x'-x|$, $\Delta y=|y'-y|$, on a une erreur initiale
$$
\eqalign{
|x'y'-xy|=|x(y'-y)+(x'-x)y'|&\leq|x|\Delta y+\Delta x|y'|\cr
&\leq|x|\Delta y+\Delta x|y|+\Delta x\Delta y.\cr
}$$
A cette erreur s'ajoute une erreur d'arrondi
$$
\Delta(x'y')\leq\varepsilon|x'y'|\leq\varepsilon(|x|+\Delta x)(|y|+\Delta y).
$$
En n\'egligeant les termes $\Delta x\Delta y$, $\varepsilon\Delta x$,
$\varepsilon\Delta y$, on obtient  la formule approximative

\bigskip
\majdeb
\hfill$
\Delta(xy)\leq|x|\Delta y+\Delta x|y|+\varepsilon|xy|.
$\hfill($*$)
\par\majfin\par

\bigskip
\noindent Soit plus g\'en\'eralement des r\'eels $x_1,\dots,x_k$, suppos\'es
repr\'esent\'es sans erreur. La formule ($*$) entra\^{\i}ne
$$
\Delta(x_1x_2\dots x_k)\leq\Delta(x_1\dots x_{k-1})|x_k|+\varepsilon|x_1\dots
x_{k-1}\cdot x_k|,
$$
d'o\`u par une r\'ecurrence ais\'ee~:
$$
\Delta(x_1x_2\dots x_k)\leq (k-1)\varepsilon|x_1x_2\dots x_k|.
$$
L'erreur sur un quotient est donn\'ee de m\^eme par $\Delta
(x/y)\leq\varepsilon|x/y|$.
On en d\'eduit pour tous exposants $\alpha_i\in\bZ$ la formule g\'en\'erale
$$
\Delta(x_1^{\alpha_1}x_2^{\alpha_2}\dots
x_k^{\alpha_k})\leq(|\alpha_1|+\dots+|\alpha_k|-1)\varepsilon|x_1^{\alpha_1}
x_2^{\alpha_2}\dots x_k^{\alpha_k}|~;
$$
on observera que $|\alpha_1|+\dots+|\alpha_k|-1$ est exactement le nombre
d'op\'erations
requises pour calculer $x_1^{\alpha_1}x_2^{\alpha_2}\dots x_k^{\alpha_k}$ par
multiplications ou divisions successives des $x_i$.

Contrairement au cas de l'addition, la majoration de l'erreur d'un  produit
{\it ne d\'epend pas de l'ordre des facteurs.}

\subsection 1.5. R\`egle de H\"orner|

\noindent On s'int\'eresse ici au probl\`eme de l'\'evaluation d'un polyn\^ome
$$
P(x)=\sum_{k=0}^na_kx^k.
$$
La m\'ethode la plus \og na\"{\i}ve\fg\ qui vient \`a l'esprit consiste \`a
poser $x^0=1$, $s_0=a_0$, puis \`a calculer par r\'ecurrence
$$
\cases{
x^k=x^{k-1}\cdot x\hfill\cr
\noalign{\vskip3pt}
u_k=a_k\cdot x^k\hfill\cr
\noalign{\vskip3pt}
s_k=s_{k-1}+u_k\hfill\cr
}\qquad{\rm pour}\quad k\geq 1.
$$
Pour chaque valeur de $k$, deux multiplications et une addition sont donc
n\'ecessaires.
Il existe en fait une m\'ethode plus efficace~:

\medskip
\state R\`egle de H\"orner -- | On factorise $P(x)$ sous la forme~:
$$
P(x)=a_0+x(a_1+x(a_2+\dots+x(a_{n-1}+xa_n)\dots)).
$$
Si l'on pose
$$
p_k=a_k+a_{k+1}x+\dots+a_nx^{n-k},
$$
cette m\'ethode revient \`a calculer $P(x)=p_0$ par r\'ecurrence descendante~:
$$
\cases{
p_n=a_n\hfill\cr
\noalign{\vskip3pt}
p_{k-1}=a_{k-1}+xp_k,\hfill&\quad$1\leq k\leq n$.\cr
}
$$
On effectue ainsi seulement une multiplication et une addition \`a chaque
\'etape, ce qui
\'economise une multiplication et donc une fraction substantielle du temps
d'ex\'ecution.

\medskip
Comparons maintenant les erreurs d'arrondi dans chacune des deux m\'ethodes, en
supposant que les r\'eels $x,a_0,a_1,\dots,a_n$ sont repr\'esent\'es sans
erreur.

\medskip
\noindent$\bullet$ {\bf M\'ethode \og na\"{\i}ve\fg.}
On a  ici $P(x)=s_n$ avec
$$
\Delta(a_k\cdot x^k)\leq
k\varepsilon|a_k||x|^k,\qquad\qquad\qquad\qquad\qquad\qquad
\qquad\qquad\quad
$$
\vskip-18pt
$$
\eqalign{
\Delta s_k&\leq \Delta
s_{k-1}+k\varepsilon|a_k||x|^k+\varepsilon(|s_{k-1}|+|u_k|)\cr
&\leq\Delta
s_{k-1}+k\varepsilon|a_k||x|^k+\varepsilon(|a_0|+|a_1||x|+\dots+|a_k||x|^k).\cr
}
$$
Comme $\Delta s_0=0$, il vient apr\`es sommation sur $k$~:
$$
\eqalign{
\Delta s_n
&\leq\sum_{k=1}^nk\varepsilon|a_k||x|^k+\varepsilon\sum_{k=1}^n
(|a_0|+|a_1||x|+\dots+|a_k||x|^k)\cr
&\leq\sum_{k=1}^nk\varepsilon|a_k||x|^k+\varepsilon\sum_{k=0}^n
(n+1-k)|a_k||x|^k.}
$$
On obtient par cons\'equent
$$
\Delta P(x)\leq (n+1)\varepsilon\sum_{k=0}^n|a_k||x|^k.
$$
$\bullet$ {\bf R\`egle de H\"orner.} Dans ce cas, on a
$$
\eqalign{
\Delta p_{k-1}&\leq\Delta(xp_k)+\varepsilon(|a_{k-1}|+|xp_k|)\cr
&\leq(|x|\Delta p_k+\varepsilon|xp_k|)+\varepsilon(|a_{k-1}|+|xp_k|)\cr
&=\varepsilon(|a_{k-1}|+2|x||p_k|)+|x|\Delta p_k.\cr
}
$$
En d\'eveloppant $\Delta P(x)=\Delta p_0$, il vient
$$
\Delta p_0\leq\varepsilon(|a_0|+2|x||p_1|)+|x|\Big(\varepsilon|a_1|+
2|x||p_2|+|x|\Big(\varepsilon|a_2|+\dots\Big)\Big)
$$
d'o\`u
$$
\eqalign{
\Delta P(x)&\leq\varepsilon\sum_{k=0}^n|a_k||x|^k+2\varepsilon
\sum^n_{k=1}|x|^k|p_k|,\cr
\Delta P(x)&\leq\varepsilon\sum_{k=0}^n|a_k||x|^k+2\varepsilon
\sum^n_{k=1}(|a_k||x|^k+\dots+|a_n||x|^n),\cr
\Delta P(x)&\leq \varepsilon\sum_{k=0}^n(2k+1)|a_k||x|^k.\cr
}
$$
On voit que la somme des coefficients d'erreur affect\'es aux termes
$|a_k||x|^k$, soit
$\displaystyle\varepsilon\sum_{k=0}^n(2k+1)=\varepsilon(n+1)^2$, est la m\^eme
que pour la
m\'ethode na\"{\i}ve~; comme $2k+1\leq 2(n+1)$, l'erreur commise sera dans le
pire des cas
\'egale \`a 2 fois celle de la m\'ethode na\"{\i}ve. N\'eanmoins, les petits
coefficients
portent sur les premiers termes calcul\'es, de sorte que la pr\'ecision de la
m\'ethode de
H\"orner sera nettement meilleure si le terme $|a_k||x|^k$ d\'ecro\^{\i}t
rapidement~:
c'est le cas par exemple si $P(x)$ est le d\'ebut d'une s\'erie convergente.

\bigskip
\majedeb
\state Exercice -- |  {\it Evaluer dans les deux cas l'erreur commise
sur les sommes partielles de la s\'erie exponentielle
$$
\sum_{k=0}^n\frac{x^k}{k!},\quad x\geq 0
$$
en tenant compte du fait qu'on a une certaine erreur d'arrondi sur
$a_k=\dfrac{1}{k!}$.
\medskip
\noindent{\rm R\'eponse.} On trouve $\Delta P(x)\leq\varepsilon(1+(n+x)e^x)$
pour la
m\'ethode na\"{\i}ve, tandis que la factorisation
$$
P(x)=1+x\Big(1+\frac{x}{2}\Big(1+\frac{x}{3}\Big(1+\dots
\Big(1+\frac{x}{n-1}\Big(1+\frac{x}{n}\Big)\Big)\dots\Big)\Big)
$$
donne $\Delta P(x)\leq\varepsilon(1+3xe^x)$, ce qui est nettement meilleur en
pratique
puisque $n$ doit \^etre choisi assez grand.}
\par\majefin\par

\subsection 1.6. Cumulation d'erreurs d'arrondi al\'eatoires|

\noindent Les majorations d'erreurs que nous avons donn\'ees plus haut
p\^echent en
g\'en\'eral par exc\`es de pessimisme, car nous n'avons tenu compte que de la
valeur
absolue des erreurs, alors qu'en pratique elles sont souvent de signe
al\'eatoire et se
compensent donc partiellement entre elles.

Supposons par exemple qu'on cherche \`a calculer une somme $s_n$ de
rang \'elev\'e d'une s\'erie convergente $S=\sum_{k=0}^{+\infty}u_k$,
les $u_k$ \'etant des r\'eels $\geq 0$ suppos\'es repr\'esent\'es sans
erreur. On pose donc
$$
\eqalign{
s_k&=s_{k-1}+u_k,\quad  s_0=u_0,\qquad\qquad\enspace\kern2pt\cr
}$$
et les erreurs $\Delta s_k$ v\'erifient
$$
\eqalign{
\Delta s_k&=\Delta s_{k-1}+\alpha_k\cr
{\rm avec}\quad \Delta s_0&=0\quad{\rm et}\quad
|\alpha_k|\leq\varepsilon(s_{k-1}+u_k)=\varepsilon s_k\leq \varepsilon S.\cr}
$$
On en d\'eduit donc
$$
\Delta s_n=\alpha_1+\alpha_2+\dots+\alpha_n
$$
et en particulier $|\Delta s_n|\leq n\varepsilon S$. Dans le pire des cas,
l'erreur est
donc proportionnelle \`a $n$. On va voir qu'on peut en fait esp\'erer beaucoup
mieux sous des hypoth\`eses raisonnables.
\bigskip\medskip

\vbox{
\state Hypoth\`eses|
\smallskip

{\parindent1.6pc
\item{(1)} Les erreurs $\alpha_k$ sont des variables al\'eatoires globalement
ind\'ependantes les unes des autres (lorsque les $u_k$ sont choisis
al\'eatoirement).

\item{(2)} L'esp\'erance math\'ematique $E(\alpha_k)$ est nulle, ce qui
signifie que les erreurs d'arrondi n'ont aucune tendance \`a se faire par 
exc\`es ou par d\'efaut.\vskip0pt
}}

\bigskip
\noindent L'hypoth\`ese (2) entra\^{\i}ne $E(\Delta s_n)=0$ tandis que
l'hypoth\`ese (1)
donne $$
{\rm var}(\Delta s_n)={\rm var}(\alpha_1)+\dots+{\rm var}(\alpha_n).
$$
Comme $E(\alpha_k)=0$ et $|\alpha_k|\leq \varepsilon S$, on a ${\rm
var}(\alpha_k)\leq
\varepsilon^2 S^2$, d'o\`u
$$
\sigma(\Delta s_n)=\sqrt{{\rm var}(\Delta s_n)}\leq \sqrt{n}\varepsilon S
$$
L'erreur quadratique moyenne cro\^{\i}t seulement dans ce cas comme $\sqrt{n}$.
D'apr\`es l'in\'egalit\'e de Bienaym\'e-Tchebychev on a~:
$$
P(|\Delta s_n|\geq \alpha\,\sigma(\Delta s_n))\leq \alpha^{-2}.
$$
La probabilit\'e que l'erreur d\'epasse $10\sqrt{n}\varepsilon S$ est donc
inf\'erieure
\`a 1\%.

\section 2. Ph\'enom\`enes de compensation|

\noindent Les ph\'enom\`enes de compensation se produisent lorsqu'on
tente d'effectuer des soustractions de valeurs tr\`es voisines. Ils
peuvent conduire \`a des pertes importantes de pr\'ecision.

\noindent Les exemples suivants illustrent les difficult\'es pouvant se
pr\'esenter et les rem\`edes \`a  apporter dans chaque cas.

\subsection 2.1. Exemple~: r\'esolution de l'\'equation $x^2-1634x+2=0$|

\noindent Supposons que les calculs soient effectu\'es avec 10 chiffres
significatifs. Les formules habituelles donnent alors
$$\eqalign{
\Delta'&=667\mil 487,\qquad\sqrt{\Delta '}\simeq 816,9987760\hfill\cr
x_1&=817+\sqrt{\Delta'}\simeq 1633,998776,\hfill\cr
x_2&=817-\sqrt{\Delta'}\simeq\hfill 0,0012240.\cr
}$$

On voit donc qu'on a une perte de 5 chiffres significatifs sur $x_2$ si l'on
effectue la
soustraction telle qu'elle se pr\'esente naturellement~!
Ici, le rem\`ede est simple~: il suffit d'observer que $x_1x_2=2$, d'o\`u
$$
x_2=\frac{2}{x_1}\simeq 1,223991125\cdot 10^{-3}.
$$
C'est donc l'{\it algorithme num\'erique utilis\'e qui doit \^etre modifi\'e.}

\subsection 2.2. Exemple~: calcul approch\'e de $e^{-10}$|

\noindent Supposons qu'on utilise pour cela la s\'erie
$$
e^{-10}\simeq\sum_{k=0}^n(-1)^k\frac{10^k}{k!},
$$
les calculs \'etant toujours effectu\'es avec 10 chiffres significatifs. Le
terme
g\'en\'eral $|u_k|={10^k}/{k!}$ est tel que
$$
\frac{|u_k|}{|u_{k-1}|}=\frac{10}{k}\geq 1\quad\hbox{d\`es que}\quad k\leq 10.
$$
On a donc 2 termes de valeur absolue maximale
$$
|u_9|=|u_{10}|=\frac{10^{10}}{10!}\simeq 2,755\cdot 10^{3}
$$
tandis que $e^{-10}\simeq 4,5\cdot 10^{-5}$. Comparons $u_{10}$ et $e^{-10}$~:

$$
\offinterlineskip\halign{
#\hfill\quad&
\tv\kern6pt\hfill#\hfill\kern6pt&
\tv\kern6pt\hfill#\hfill\kern6pt&
\tv\kern6pt\hfill#\hfill\kern6pt&
\tv\kern6pt\hfill#\hfill\kern6pt&
\tv\kern6pt\hfill#\hfill\kern6pt&
\tv\kern6pt\hfill#\hfill\kern6pt&
\tv\kern6pt\hfill#\hfill\kern6pt&
\tv\kern6pt\hfill#\hfill\kern6pt&
\tv\kern6pt\hfill#\hfill\kern6pt&
\tv\kern6pt\hfill#\hfill\kern6pt&
\tv#\cr
&\multispan{10}\hrulefill&\omit\cr
$u_{10}$~:&2&7&5&5{\bf ,}&$\cdot$&$\cdot$&$\cdot$&$\cdot$&$\cdot$&$\cdot$&\cr
&\multispan{10}\hrulefill&\omit\cr
%&\multispan4\hfill ,&\omit&\omit&\omit&\omit&\omit&\omit&\omit\cr
\noalign{\vskip6pt}
&\multispan{10}\hrulefill&\omit\cr
$e^{-10}~:$&&&&0{\bf ,}&0&0&0&0&4&5&\cr
&\multispan{10}\hrulefill&\omit\cr
%&\multispan4\hfill ,&\omit&\omit&\omit&\omit&\omit&\omit&\omit\cr
}
$$

\noindent Ceci signifie qu'au moins 8 chiffres significatifs vont \^etre perdus
par
compensation des termes $u_k$ de signes oppos\'es. Un rem\`ede simple consiste
\`a
utiliser la relation
$$
e^{-10}=1/e^{10}\quad{\rm avec}\quad e^{10}\simeq\sum_{k=0}^n\frac{10^k}{k!}.
$$
On essaiera dans la mesure du possible d'{\it \'eviter les sommations dans
lesquelles des termes de signes oppos\'es se compensent.}

\subsection 2.3. Exemple~: calcul approch\'e de $\pi$ par les
polyg\^ones inscrits|

\noindent Soit $P_n$ le demi-p\'erim\`etre du polygone r\'egulier \`a
$n$ c\^ot\'es inscrit dans un cercle de rayon 1.

\noindent Le c\^ot\'e de ce polygone vaut $2\cdot
R\,\sin\,\pi/n=2\,\sin\,\pi/n$, d'o\`u
$$
\eqalign{
P_n&=n\,\sin\frac{\pi}{n},\cr
P_n&=\pi-\frac{\pi^3}{6n^2}+o\Big(\frac{1}{n^3}\Big).\cr
}
$$
On obtient donc une approximation de $\pi$ avec une pr\'ecision de
l'ordre de $6/n^2$.

\InsertFig 0.000 46.000 {
1.2 mm unit  /nbpt 10 def
-35.000   0.000 translate
 70.000  20.000 moveto
[  70.000  20.000  110.000  10.000  110.000  30.000   70.000  20.000
] polygon
stroke
108.492   5.224 moveto
 70.000  20.000  41.231 -21.000  21.000 arc stroke
 80.000  20.000 moveto
 70.000  20.000  12.500   0.000 -14.306 arcn stroke
 80.700  20.000 moveto
 70.000  20.000  13.200   0.000 -14.162 arcn stroke
 83.000  20.000 moveto
 70.000  20.000  15.500   0.000  14.303 arc stroke
 83.700  20.000 moveto
 70.000  20.000  16.200   0.000  14.331 arc stroke
 96.000  28.000 moveto
 14.300 rotmem
  1.000 nbpt mul it  (R = ) show currentpoint
clear
  1.000 nbpt mul rm  (1) show currentpoint
clear
rotrev
[   2.000   1.000 ] 0 setdash
 70.000  20.000 moveto
110.000  20.000 lineto
stroke}
\LabelTeX  62.000  26.000  $\pi/n$\ELTX
\EndFig

\noindent Pour \'evaluer $P_n$, on utilise une m\'ethode de dichotomie
permettant de calculer par r\'ecurrence
$$
x_k=P_{2^k}=2^k\,\sin\frac{\pi}{2^k}.
$$
Si $\alpha$ est un angle compris entre 0 et $\frac{\pi}{2}$ on a
$$
\sin\,\frac{\alpha}{2}=\sqrt{\frac{1}{2}(1-\cos\,\alpha)}=
\sqrt{\frac{1}{2}(1-\sqrt{1-\sin^2\alpha)}}.\eqno(*)
$$
En substituant $\alpha={\pi}/{2^k}$, on en d\'eduit les formules
$$
\cases{
x_{k+1}=2^k\sqrt{2(1-\sqrt{1-(x_k/2^k)^2)}}\hfill\cr
\noalign{\vskip3pt}
x_1=2,\hfill\cr
}
$$
et d'apr\`es ce qu'on a dit plus haut
$\displaystyle\lim_{k\rightarrow+\infty}x_k=\pi$.

\noindent Ce n'est pourtant pas du tout ce qu'on va observer sur machine~! 
D\`es que $(x_k/2^k)^2$ sera inf\'erieur \`a la pr\'ecision relative des
calculs, l'ordinateur va
donner
$$
\sqrt{1-(x_k/2^k)^2}=1\quad\hbox{d'o\`u}\quad x_{k+1}=0.
$$
Pour \'eviter cette difficult\'e, il suffit de remplacer $(*)$ par
$$
\sin\,\frac{\alpha}{2}=\sqrt{\frac{1}{2}\
\frac{1-\cos^2\alpha}{1+\cos\,\alpha}}=\frac{\sin\,\alpha}
{\sqrt{2(1+\cos\,\alpha)}}
\eqno(**)
$$
d'o\`u
$$
\sin\,\frac{\alpha}{2}=\frac{\sin\,\alpha}{\sqrt{2(1+\sqrt{1-\sin^2\alpha)}}}.
$$
On obtient alors la formule de r\'ecurrence
$$
x_{k+1}=\frac{2x_k}{\sqrt{2(1+\sqrt{1-(x_k/2^k)^2})}}
$$
qui \'evite le ph\'enom\`ene de compensation pr\'ec\'edent, de sorte que le
calcul des
$x_k$ peut \^etre pouss\'e beaucoup plus loin.

On obtiendra une m\'ethode plus efficace encore en observant qu'on peut
\'evaluer
$\cos\,\alpha$ dans $(**)$ par la formule
$\cos\,\alpha=\dfrac{\sin\,2\alpha}{2\,\sin\,\alpha}$. Ceci donne
$$
\eqalign{
&\sin\,\frac{\alpha}{2}=\sin\,\alpha\sqrt{\frac{\sin\,\alpha}{2\,\sin\,\alpha
+\sin\,2\alpha}},\quad\hbox{d'o\`u}\cr
&x_{k+1}=x_k\sqrt{\frac{2x_k}{x_k+x_{k-1}}}.
}$$
Deux valeurs d'initialisation sont alors requises pour d\'emarrer, par exemple
$x_1=2$ et
$x_2=2\sqrt{2}$.

\section 3. Ph\'enom\`enes d'instabilit\'e num\'erique|

\noindent Il s'agit de ph\'enom\`enes d'amplification des erreurs d'arrondi.
Une telle
amplification se produit assez fr\'equemment dans le cas de calculs
r\'ecurrents ou
it\'eratifs.

\subsection 3.1. Cas d'un calcul r\'ecurrent|

\noindent Supposons \`a titre d'exemple qu'on cherche \`a \'evaluer
num\'eriquement l'int\'egrale
$$
I_n=\int^1_0\frac{x^n}{10+x},\qquad n\in\bN.
$$
Un calcul imm\'ediat donne
$$
\eqalign{
I_0&=\int^1_0\frac{dx}{10+x}=\Big[\ln\,(10+x)\Big]^1_0=\ln\,\frac{11}{10},\cr
I_n&=\int^1_0\frac{x}{10+x}\cdot
x^{n-1}dx=\int^1_0\Big(1-\frac{10}{10+x}\Big)x^{n-1}dx\cr
&=\int^1_0x^{n-1}dx-10\int^1_0\frac{x^{n-1}}{10+x}\,dx
=\frac{1}{n}-10\,I_{n-1}.\cr
}
$$
Ceci permet de calculer $I_n$ par r\'ecurrence avec
$$
\cases{
I_0=\ln\,\dfrac{11}{10}\hfill\cr
\noalign{\vskip3pt}
I_n=\frac{1}{n}-10\,I_{n-1}.\hfill\cr
}
$$
Ce probl\`eme apparemment bien pos\'e math\'ematiquement conduit
num\'eriquement \`a
des r\'esultats catastrophiques. On a en effet
$$
\Delta I_n\simeq10\,\Delta I_{n-1},
$$
m\^eme si on n\'eglige l'erreur d'arrondi sur $1/n$. L'erreur sur $I_n$ explose
donc
exponentiellement, l'erreur initiale sur $I_0$ \'etant multipli\'ee par $10^n$
\`a
l'\'etape $n$. Comment faire alors pour calculer par exemple $I_{36}$~?
La suite $x^n$ \'etant d\'ecroissante pour $x\in[0,1]$, on voit que la suite
$I_n$ est elle-m\^eme d\'ecroissante. Comme $10\leq  10+x\leq 11$, on a de plus
$$
\frac{1}{11(n+1)}\leq I_n\leq\frac{1}{10(n+1)}.
$$
L'approximation $I_n\simeq \dfrac{1}{11(n+1)}$ donne une erreur absolue
$\leq\dfrac{1}{110(n+1)}$ et donc une erreur relative $\dfrac{\Delta
I_n}{I_n}\leq\dfrac{1}{10}$. Ceci donne l'ordre de grandeur mais n'est pas
tr\`es
satisfaisant. L'id\'ee est alors de {\it renverser la r\'ecurrence} en posant
$$
I_{n-1}=\frac{1}{10}\Big(\frac{1}{n}-I_n\Big).
$$
En n\'egligeant l'erreur sur $\dfrac{1}{n}$, on a donc cette fois $\Delta
I_{n-1}\simeq\dfrac{1}{10}\,\Delta I_n$, estimation qui va dans le bon sens.
Si l'on part de $I_{46}\simeq\dfrac{1}{11.47}$, on obtiendra pour
$I_{36}$ une erreur relative sans doute meilleure que $10^{-10}$.

\bigskip
\majedeb
\state Exercice -- |  {\it Montrer que $0\leq
I_n-I_{n+1}\leq\dfrac{1}{10(n+1)(n+2)}$, et en d\'eduire \`a partir de
la formule exprimant $I_n$ en fonction de $I_{n+1}$ que l'on a en fait
l'estimation
$$
\frac{1}{11(n+1)}\leq I_n\leq\frac{1}{11(n+1)}+\frac{1}{110(n+1)(n+2)},
$$
donc $\Delta I_n\simeq\dfrac{1}{11(n+1)}$ avec erreur relative $\leq
\dfrac{1}{10(n+2)}$.}
\par\majefin\par

\bigskip
On voit donc le r\^ole fondamental jou\'e par le coefficient d'amplification de
l'erreur, 10 dans le premier cas, $1/10$ dans le second. En g\'en\'eral, si on
a un
coefficient d'amplification $A>1$, il est imp\'eratif de limiter le nombre $n$
d'\'etapes en sorte que $A^n\varepsilon$ reste tr\`es inf\'erieur \`a 1, si
$\varepsilon$ est la pr\'ecision relative des calculs.

\subsection 3.2. Calculs it\'eratifs|

\noindent Soit \`a calculer une suite $(u_n)$ d\'efinie par sa valeur
initiale $u_0$ et par la relation de r\'ecurrence
$$
u_{n+1}=f(u_n),
$$
o\`u $f$ est une fonction donn\'ee. On a donc $u_n=f^n(u_0)$ o\`u $f^n=f\circ
f\circ\dots \circ f$ est la $n$-i\`eme it\'er\'ee de $f$. On consid\`ere par
exemple la
suite $(u_n)$ telle que
$$
u_0=2,\qquad u_{n+1}=|\ln\,(u_n)|,
$$
dont on cherche \`a \'evaluer le terme $u_{30}$. Un calcul effectu\'e en
Turbo-Pascal
donne imm\'ediatement $u_{30}\simeq 0,880833175$.

A la lumi\`ere de l'exemple pr\'ec\'edent, il est n\'eanmoins l\'egitime de se
demander
si ce calcul est bien significatif, compte tenu de la pr\'esence des erreurs
d'arrondi.
En partant de valeurs de $u_0$ tr\`es voisines de 2, on obtient en fait les
r\'esultats
suivants (arrondis \`a $10^{-9}$ pr\`es)~:
$$
\vbox{\offinterlineskip\halign{
\hfill#\hfill\quad&
\tv\quad\hfill#\hfill\quad&
\tv\quad\hfill#\hfill\quad&
\tv\quad\hfill#\quad&
\tv\kern2pt\tv\quad\hfill#\cr
$u_0$&2,000000000&2,000000001&1,999999999&$5\cdot 10^{-10}$\cr
$u_5$&5,595485181&5,595484655&5,595485710&$9\cdot 10^{-8}$\cr
$u_{10}$&0,703934587&0,703934920&0,703934252&$5\cdot 10^{-7}$\cr
$u_{15}$&1,126698502&1,126689382&1,126707697&$8\cdot 10^{-6}$\cr
$u_{20}$&1,266106839&1,266256924&1,265955552&$10^{-4}$\cr
$u_{24}$&1,000976376&1,001923276&1,000022532&$10^{-3}$\cr
$u_{25}$&0,000975900&0,001921429&0,000022532&100\%\cr
$u_{26}$&6,932150628&6,254686211&10,700574400&50\%\cr
$u_{30}$&0,880833175&0,691841353&1,915129896&100\%\cr
}}
$$
La derni\`ere colonne donne l'ordre de grandeur de l'\'ecart relatif $\Delta
u_n/u_n$
observ\'e entre la deuxi\`eme ou troisi\`eme colonne et la premi\`ere colonne.
On voit
que cet \'ecart augmente constamment pour atteindre environ $10^{-3}$ sur
$u_{24}$. Pour
le calcul de $u_{25}$, il se produit une v\'eritable catastrophe num\'erique~:
l'\'ecart
relatif devient voisin de 100\%~! Il en r\'esulte que toutes les valeurs
calcul\'ees \`a
partir de $u_{25}$ sont certainement non significatives pour une pr\'ecision
des calculs
de $10^{-9}$.

Pour comprendre ce ph\'enom\`ene, il suffit d'observer qu'une erreur $\Delta x$
sur la
variable $x$ entra\^{\i}ne une erreur $\Delta f(x)$ sur $f(x)$,
approximativement donn\'ee
par

\bigskip
\majdeb
\centerline{$
\Delta f(x)=|f'(x)|\,\Delta x.
$}
\par\majfin\par

\bigskip
Ceci se voit bien s\^ur en approximant $f(x+\Delta x)-f(x)$ par sa
diff\'erentielle
$f'(x)\Delta x$, lorsque $f$ est d\'erivable au point $x$. Le coefficient
d'amplification
de l'erreur absolue est donc donn\'e par la valeur absolue de la d\'eriv\'ee
$|f'(x)|$~;
ce coefficient peut \^etre parfois assez grand. Souvent dans les calculs
num\'eriques (et
ici en particulier), il est plus pertinent de consid\'erer les erreurs
relatives. La formule
$$
\frac{\Delta f(x)}{|f(x)|}=\frac{|f'(x)||x|}{|f(x)|}\ \frac{\Delta x}{|x|}
$$
montre que le coefficient d'amplification de l'erreur relative est
$|f'(x)||x|/|f(x)|$.
Dans le cas $f(x)=\ln\,(x)$ qui nous int\'eresse, ce coefficient vaut
$1/|\ln\,x|$~; il
devient tr\`es grand lorsque $x$ est proche de 1, comme c'est le cas par
exemple pour
$u_{24}$.

\section 4. Probl\`emes|

{\parindent1.4pc

\noindent{\bf 4.1.} Soit $x\geq 0$~; on note
$\displaystyle F(x)=\dfrac{2}{\sqrt{\pi}}\int^x_0e^{-t^2}dt$.
\medskip

\item{(a)} Encadrer $F(x)$ par deux entiers cons\'ecutifs.
\medskip

\item{(b)} En rempla\c cant $e^{-t^2}$ par un d\'eveloppement en s\'erie
enti\`ere de
$x$, exprimer $F(x)$ comme somme d'une s\'erie. On choisit $x=3$~; calculer les
10
premiers termes de la s\'erie. En d\'eduire que pour $x\geq 3$ on a un
ph\'enom\`ene de
compensation dans le calcul de la somme des premiers termes de la s\'erie.

\medskip
\item{(c)} On d\'efinit $g(x)$ par $F(x)=e^{-x^2}g(x)$.

Montrer que $g$ est solution d'une \'equation diff\'erentielle.

Exprimer $g(x)$ comme somme d'une s\'erie enti\`ere en $x$.

\medskip
\item{(d)} En d\'eduire l'expression $F(x)=\sum^{+\infty}_{n=0}a_n(x)$
o\`u les $a_n(x)$ sont tous positifs. D\'eterminer $a_0(x)$ et donner la
solution de
r\'ecurrence entre $a_n(x)$ et $a_{n-1}(x)$. Montrer l'in\'egalit\'e
$$
\sum_{n=N+1}^{+\infty}a_n(x)\leq a_N\ \frac{x^2}{N-x^2}\quad ({\rm pour}\
N>x^2)
$$
\item{(e)} En utilisant les r\'esultats pr\'ec\'edents, \'ecrire un programme
en langage
Pascal qui, \`a la lecture de $x$ et d'un entier $k\leq 1$ calcule une valeur
approch\'ee de $F(x)$ \`a $10^{-k}$ pr\`es.


\bigskip
\noindent{\bf 4.2.} Soit $(I_n)_{n\in\bN}$ la suite des int\'egrales
$$
I_n=\int^1_0\frac{x^n}{6+x-x^2}\ dx.
$$
\item{(a)} Montrer que $I_n$ v\'erifie une relation de r\'ecurrence de la forme
$$
I_{n+1}=\alpha I_n+\beta I_{n-1}+c_n\eqno(*)
$$
\item{} o\`u $\alpha,\beta$ sont des constantes et $(c_n)$ une suite
num\'erique
explicite.

\medskip
\item{(b)} On envisage le calcul r\'ecurrent de $I_n$ \`a partir de $I_0$ et
$I_1$ par
la formule $(*)$.

\item{} On suppose que les valeurs de $I_0$ et $I_1$ sont affect\'ees d'erreurs
d'arrondis
$\varepsilon_0$ et $\varepsilon_1$, et on note $\varepsilon_n$ l'erreur qui en
r\'esulte
sur $I_n$ (on n\'eglige ici l'erreur sur le calcul de $c_n$ et les erreurs
d'arrondi
pouvant intervenir dans l'application de la formule $(*)$).

\smallskip
\itemitem{$(\alpha)$} D\'eterminer 
$\varepsilon_n$ en fonction de $\varepsilon_0$ et $\varepsilon_1$.

\smallskip
\itemitem{$(\beta)$} Est-il possible de calculer $I_{50}$ par ce proc\'ed\'e
avec un ordinateur donnant une pr\'ecision relative de $10^{-10}$~?


\bigskip
\noindent{\bf 4.3.} Etant donn\'e une suite $x_k$, $k=1,\dots,n$ de r\'eels, 
on note $\mu_n=\dfrac{1}{n}\ \sum^n_{k=1}x_k$\break la moyenne
et $\sigma_n=\sqrt{\sigma^2_n}$ l'\'ecart type avec
$\sigma^2_n=\dfrac{1}{n}\ \sum^n_{k=1}(x_k-\mu_n)^2$.
\medskip

\item{(a)} Soit  $q_n=\sum^n_{k=1}x^2_k$. Exprimer $\sigma^2_n$ en fonction
de $q_n$ et de $\mu_n$.

\medskip
\item{(b)} Ecrire un programme qui calcule les moyennes et l'\'ecart type d'un
nombre
ind\'etermin\'e de r\'eels. Les donn\'ees r\'eelles sont entr\'ees au clavier~;
apr\`es
chaque entr\'ee on affichera la moyenne et l'\'ecart type de la suite des
nombres
d\'ej\`a entr\'es.

\medskip
\item{(c)} On suppose que pour $k=1,\dots,n$ on a $x_k=\mu+\varepsilon_k$ avec
$|\varepsilon_k|<\varepsilon$ o\`u $\varepsilon$ est petit devant $\mu$.
Montrer que l'on a l'in\'egalit\'e $\big|\dfrac{q_n}{n}-\mu^2\big|\leq 
3\mu\varepsilon$.

\item{} En d\'eduire que la m\'ethode de calcul de $\sigma_n$ utilisant la
formule du (a)
est inadapt\'ee pour une telle suite.

\medskip

\item{(d)} On veut obtenir un algorithme de calcul de $\sigma_n$ plus stable.

\item{} Etablir les \'egalit\'es~:
$$
\eqalign{
&(n+1)\sigma^2_{n+1}=n\sigma^2_n+n(\mu_{n+1}-\mu_n)^2+(x_{n+1}-\mu_{n+1})^2,\cr
&(n+1)\mu_{n+1}=n\mu_n+x_{n+1}.\cr
}
$$
\item{} En d\'eduire $\sigma^2_{n+1}=\dfrac{n}{n+1}\ \sigma^2_n+\dfrac{1}{n}\
(x_{n+1}-\mu_{n+1})^2$.

\medskip
\item{(e)} Reprendre la question (b) avec le nouvel algorithme.

\medskip
\item{(f)} On consid\`ere une suite de r\'eels $x_k=1+\varepsilon\
\dfrac{2k-n-1}{n-1}$, $k=1,\dots,n$. D\'eterminer sa moyenne et son
\'ecart type.

\medskip \item{(g)} M\^eme question pour la suite des $2^n$ r\'eels
$x_k$, $k=1,\dots,2^n$ telle que pour $p=0,\dots,n$ on ait $\rmC^p_n$
termes \'egaux \`a $\mu+\dfrac{2p-n}{\sqrt{n}}$ (On pourra remarquer
que $p^2\rmC^p_n=pn\rmC^{p-1}_{n-1}$).
\par}

\bye

