\documentclass[11pt,a4paper]{article}
\parindent=0cm
 
   \textheight=23,5cm
   \textwidth=15,7cm
\oddsidemargin=0cm

\usepackage[frenchb]{babel}
\usepackage[utf8]{inputenc}


\usepackage{url,amsmath,amssymb,amsthm,amsfonts,epsfig,eurosym,geometry}
\selectlanguage{french}

\def\var{\mathrm{var}}
\def\BB{\mathcal{B}}
\def\FF{\mathcal{F}}
\def\Pa{\mathcal{P}}
\def\M{\mathcal{M}}
\def\P{\mathcal{P}}
\def\T{\mathcal{T}}
\def\L{\mathcal{L}}

\def\AA{\mathfrak{A}}
\def\S{\mathfrak{S}}
\def\SS{\mathcal{S}}
\def\CC{\mathcal{C}}

\def\A{\mathbb{A}}
\def\B{\mathbb{B}}
\def\D{\mathbb{D}}
\def\U{\mathbb{U}}


\def\bq{\begin{quote}\begin{em}}
\def\eq{\end{em}\end{quote}}

\def\un{\mathbf{1}}
\def\dd{\mathrm{d}}
\def\ee{\mathrm{e}}
\def\ii{\mathrm{i}}
\def\eps{\varepsilon}

\def\hth{\widehat{\theta}}
\def\ovx{\overline{X}}
\def\ima{\mathrm{Im}}

\newcommand{\C}{\mathbb C}
\newcommand{\R}{\mathbb R}
\newcommand{\N}{\mathbb N}
\newcommand{\Z}{\mathbb Z}
\newcommand{\Q}{\mathbb Q}

\newcommand{\espmes }{(\Omega, {\cal F})}
\newcommand{\espprob}{(\Omega, {\cal F}, P)}
\newcommand{\va}{variable al\'eatoire }
\newcommand{\vas}{variables al\'eatoires }
\newcommand{\Vas}{Variables al\'eatoires }
\newcommand{\cov}{\hbox{cov}}
\newtheorem{theo}{Théorème}[section]
\newtheorem{coro}{Corollaire}[section]
\newtheorem{lemm}{Lemme}[section]
\newtheorem{defi}[theo]{Définition}
\newtheorem{prop}[theo]{Proposition}
\newtheorem{rque}[theo]{Remarque}
\newtheorem{exem}[theo]{Exemple}
\newtheorem{exer}[theo]{Exercice}
\def\card{\mathrm{card}\,}

\newcommand{\tran}[1]{\,{}^{t}\!#1}
\renewcommand\det{\mathop{\textrm{d\'et}}}

\renewcommand{\le}{\leqslant}
\renewcommand{\ge}{\geqslant}
\def\demi{\frac{1}{2}}
\def\ds{\displaystyle}
\def\btheo{\begin{theo}}\def\etheo{\end{theo}}
\def\bcoro{\begin{coro}}\def\ecoro{\end{coro}}
\def\blemm{\begin{lemm}}\def\elemm{\end{lemm}}
\def\bdefi{\begin{defi}}\def\edefi{\end{defi}}
\def\bprop{\begin{prop}}\def\eprop{\end{prop}}
\def\brque{\begin{rque}}\def\erque{\end{rque}}
\def\bexem{\begin{exem}}\def\eexem{\end{exem}}

\def\bexer{\begin{exer}\begin{em}}\def\eexer{\end{em}\end{exer}}
\usepackage{graphics}
\usepackage{color}              % Need the color package
\usepackage{epsfig}
\def\F{{\cal F}}
\def\1{{\rm 1\kern-.8ex 1}}
\def\euro{\mbox{\raisebox{.25ex}{{\it =}}\hspace{-.5em}{\sf C}}} 
\newcommand{\bin}[2]{\left( \begin{array}{c} 
#1 \\#2\end{array}\right)}
\begin{document}

\noindent Intro mod\'elisation num\'erique \hfill Feuille 1: Nombres flottants
\hfill {L3 2016/2017}\\


{\bf Notations}~: un nombre flottant normalis\'e $f$ est de la forme
\[ f=\pm \left(\sum_{k=1}^r\frac{a_k}{b^k} \right) b^e=\pm m b^{e-r},
\quad
m=\sum_{k=1}^r a_k b^{r-k}, \]  
\begin{itemize}
\item  $b>1$ est un entier (la base), la plupart du temps 2, le m\^eme
pour tous les flottants
\item $r$ est le nombre de chiffres (ou bits) significatifs, le m\^eme
pour tous les flottants (sauf en multi-pr\'ecision)
\item $m<b^r$ est la mantisse de $f$, un entier positif
dont l'\'ecriture en base $b$ est
$a_1,\ldots ,a_r $ avec $a_1 \in ]0,b-1]$.
(si $b=2$ on ne stocke pas $a_1$ qui est toujours \'egal '`a 1).
\item $e$ est l'exposant de $f$, $L\leq e\leq U$.
\end{itemize}
La norme IEEE-754 (langages de programmation, logiciels de calcul
num\'erique...) utilise $b=2$, $r=53$, $L=-1021$, $U=1024)$. Xcas utilise
cette norme mais avec troncature des 5 bits les moins
significatifs de mantisse \`a 0 en pr\'ecision normale ($r=48$), en
multipr\'ecision $r$ est fix\'e par l'utilisateur (\verb|Digits| ou
second argument de \verb|evalf|), et $L,l$ valent au moins $2^{30}$
en valeur absolue.

La saisie d'un nombre flottant se fait en \'ecrivant la
mantisse en base 10 (avec \'eventuellement un point d\'ecimal), suivi
\'eventuellement de la lettre \verb|e| et de l'exposant en base 10, 
m\^eme si le
flottant est stock\'e en base 2. Par exemple \verb|6.02e23| pour le
nombre d'Avogadro. Attention, en Xcas, \verb|10^10| est un entier
exact, \`a ne pas confondre avec le flottant \verb|1e10|.


\paragraph*{Exercice 1 Ecriture en base $b$ d'un rationnel}
Soit $b$ un entier supérieur ou égal à 2.
Soit $p$ et $q$ deux entiers strictement positifs tels que $p<q$.

Soit $\frac{p}{q}=\ds\sum_{k=1}^{+\infty}\frac{a_k}{b^k}$ le
développement en base $b$ de $\frac{p}{q}$.
\begin{enumerate}
 \item Montrer que la suite $(a_k)$ est obtenue de la façon suivante :
\[r_0=p,\ \ \ \ a_1\ \mbox{est le quotient de la division euclidienne de }\  r_0b\  \mbox{par }\ q\  \mbox{et }\ r_1 \ \mbox{est le reste}\]
Et si $(r_0,\cdots,r_n)$ et $(a_1,\cdots,a_n)$ sont définis alors
$a_{n+1}$ est le quotient de la 
division euclidienne de $r_nb$ par $q$ et $r_{n+1}$ est son reste.
\item Montrer que la suite $(a_k)$ est périodique à partir d'un certain rang.
\item Donner le developpement de $x=\frac{1}{10}$ en base 2. 
Soit $\widehat{x}$ le flottant arrondi
  normalis\'e avec 12=11+1 bits de mantisse correspondant \`a $x$. 
Calculer $\widehat{x}$. M\^eme question pour $y=2/5$.
Faire la somme $\widehat{x}+\widehat{y}$.
\item Expliquez le r\'esultat de l'op\'eration \verb|0.5-(0.1+0.4)| dans Xcas.
\end{enumerate}

\paragraph*{Exercice 2} Durant la premi\`ere guerre du Golfe, une batterie anti-missile Patriot a rat\'e l'interception d'un missile Scud qui causa la mort de 28 personnes. Nous allons essayer de comprendre pourquoi.\\
L'ordinateur de la batterie Patriot repr\'esente les nombres par des nombres \`a virgule fixe en base 2 avec une pr\'ecision de 23 chiffres apr\`es la virgule. L'horloge de la batterie compte le temps en dixi\`eme de seconde. Pour obtenir le temps en seconde, le programme divise le temps donn\'e par l'horloge par dix.\\
1) En utilisant l'écriture en base 2 de $\frac{1}{10}$ trouvée dans l'exercice précédent, donner une estimation de l'erreur d'approximation de 1/10 par la repr\'esentation des nombres de la batterie Patriot.\\
2) La batterie fonctionne pendant 100 heures. Donner une estimation de l'erreur introduite dans le temps en secondes.\\
3) Sachant qu'un missile Scud voyage \`a 1676 m\`etres par seconde, donner l'erreur de distance commise par la batterie anti-missile.

\paragraph*{Exercice 3}
Donner la valeur exacte de chaque expression suivante.
$$\sum_{n=1}^{+\infty} \frac 1n~~,~~~\left(10^{20}+1\right)-10^{20}~~,~~~\binom{1000}{2}=\frac{1000!}{2(998!)}~.$$
On donne ces calculs \`a faire \`a une machine. 
Sa r\'eponse correspondra-t-elle \`a l'expression exacte ? (Discuter
selon le mode de calcul, exact ou approch\'e, et selon la pr\'ecision).


\paragraph*{Exercice 4}
Les nombres suivants peuvent se calculer selon deux
mani\`eres not\'ees $A$ et $B$. Quelle est l'expression la plus judicieuse
pour le calcul approch\'e sur machine~?
\begin{itemize}
\item $A=1-\cos^2(10^{-8}), \quad B=\sin^2(10^{-8})$
\item
$A=~\sum_{n=1}^{10^{9}} \frac 1n, \quad B=\sum_{n=0}^{10^{9}-1} \frac
1{10^{9}-n}~.$
\item
$A=~e^{-10}\approx \sum_{n=0}^{30}\frac {(-10)^n}{n!},
\quad B=~\frac 1{e^{10}} \approx \left(\sum_{n=0}^{30}\frac
  {10^n}{n!}\right)^{-1}$
\end{itemize}
Proposez une m\'ethode de calcul approch\'e pr\'ecise pour calculer~:
$$ \frac{1}{\sqrt{10^{10}+1}-\sqrt{10^{10}-1}}$$

\paragraph*{Exercice 5}
On effectue la somme de $n$=100 nombres flottants de taille comparable, on suppose que chaque
op\'eration d'addition introduit une erreur absolue distribu\'ee selon
une loi uniforme sur $[-1,1]$. Effectuer plusieurs simulations de la valeur de la
somme des erreurs, observer la distribution des erreurs obtenues.
Commandes Xcas utiles: \verb|M:=ranm(n,N,uniformd,-1,1)| pour
g\'en\'erer une matrice $n$ par $N$ d'erreurs selon la loi uniforme,
\verb|sum(M)| effectue la somme de ces erreurs par colonnes et \verb|histogram|
permet de les repr\'esenter. On pourra comparer l'histogramme de
\verb|sum(M)/sqrt(n)| et la loi normale de moyenne nulle et d'\'ecart-type
\verb|stddev(uniformd,-1,1)| (=$\frac{1}{\sqrt{3}}$)
 celle de la loi uniforme consid\'er\'ee.

\paragraph*{Exercice 6 : \'evaluation d'un polyn\^ome en un r\'eel.}
On veut calculer la valeur en $x$ d'un polynome $\ds\sum_{k=0}^na_kX^k$.
On suppose que $a_0,\ldots,a_n,x$ sont repr\'esent\'es exactement
par des nombres flottants (on pourra aussi r\'efl\'echir aux erreurs
introduites si la repr\'esentation n'est pas exacte!).


\begin{enumerate}
 \item On fait le calcul de la manière suivante :
$s_0=a_0, u_1=x, v_1=a_1u_1$ et pour $k \geq 1$
$$ s_k= s_{k-1} + v_k,\ \ \ \ u_{k+1}=u_k x,\ \  \ \ v_{k+1}=
u_{k+1} a_{k+1}$$
Combien a t-on fait de multiplications et d'additions ?

\item Méthode de Hörner\\
Le principe est d'écrire 
$P(X)=\ds\sum_{k=0}^na_kX^k=a_0+X(a_1+X(\cdots +X(a_{n-1}+Xa_n)))$.
On pose pour $0\leq k\leq n$,
\[p_k=a_k+a_{k+1}x+\cdots+a_nx^{n-k}\]
On a alors $P(x)=p_0$. On calcule les termes par une récurrence descendante \\
$p_n=a_n$ et $p_{k-1}=a_{k-1}+xp_k$\\
Combien a t-on fait de multiplications et d'additions ?

\item Ecrire deux programmes évaluant un polynome donné par ses coefficients en $x$.
Comparer en utilisant \verb|time()| la rapidité d'execution.\\
Dans Xcas, la fonction \verb|horner(P,x)| évalue en x le polynome P par la méthode de Hörner.
On peut définir un polynome par son expression symbolique ou la liste
de ses coefficients par ordre d\'ecroissant. On peut faire le calcul
avec des rationnels et arrondir la valeur exacte, que gagne-t-on alors~? 
Qu'observe-t-on en terme de temps d'ex\'ecution~?

\item {\em Difficile}~: Estimation des erreurs par ces deux m\'ethodes
  (attention, question d\'elicate!).\\
On note $u$ l'erreur d'arrondi/troncature machine.
On rappelle que si $x$ et $y$ sont des flottants, si $\ast\in  \{+,\times\}$, alors il existe $\epsilon$ tel que $x\circledast y=(1+\epsilon)(x\ast y)$ et $|\epsilon|\leq u$.\\
En utilisant la premi\`ere m\'ethode, on a~:
$\widehat {s_0}=a_0; \widehat{u_1}=x; \widehat{v_1}=a_1\otimes \widehat{u_1}$ et pour $k\geq 1$,
\[\widehat{s_k}=\widehat{s_{k-1}}\oplus  \widehat{v_k},\ \ \ \widehat{u_{k+1}}=\widehat{u_k}\otimes x,\ \  \ \widehat{v_{k+1}}= \widehat{u_{k+1}}\otimes a_{k+1}\]
La valeur donnée par la machine est donc $\widehat{s_n}$.\\
On note $\gamma_n(u)=(1+u)^{n+1}-1$. Montrer par récurrence que
$$|\ds\sum_{k=0}^na_kx^k-\widehat{s_n}| \leq \gamma_n(u)\
\ds\sum_{k=0}^n|a_k||x|^k$$
(On pourra montrer
$| \hat{v}_n / v_n | \leq (1+u)^n$ ,
$| \hat{v}_n  -  v_n | \leq \gamma_{n-1}(u) |v_n|$, 
$| \hat{s}^n | \leq (1+u)^{n+1} sum | a_k x^k |$).
Peut-on en d\'eduire une information sur l'erreur relative du
r\'esultat renvoy\'e, en faisant \'eventuellement des hypoth\`eses sur
les signes~?\\
Les processeurs de calcul flottant utilisent souvent une repr\'esentation plus
pr\'ecise que le standard IEE-754, avec 64 bits de mantisse (+1
implicite) au lieu de 52 bits (+1 implicite) (ou des 47+1 bits de Xcas). Si on effectue les
calculs interm\'ediaires de cette mani\`ere (ou plus g\'en\'eralement
en multi-pr\'ecision), qu'y gagne-t-on~?\\
En notant $\widehat{p_k}$, le flottant obtenu par la méthode
de H\"orner, on peut montrer que 
\[|P(x)-\widehat{p_0}|\leq\ds\sum_{k=0}^n\gamma_{2k+2}(u)|a_k||x|^k\]
Cette majoration sera meilleure que celle obtenue par la premi\`ere
méthode si les termes $a_kx^k$ décroissent vers 0, mais elle reste
du m\^eme type.


\end{enumerate}

\paragraph*{Exercice 7}
On souhaite calculer l'int\'egrale
$$
I_n=\int_0^1\frac{x^n}{10+x}dx
$$
pour $n\in\N$.
\begin{enumerate}
\item Calculer $I_0$. % Montrer que $I_0=\mbox{ln}(11/10)$.
\item Montrer que pour tout $n\in\N^*$,
$$I_n=\frac{1}{n} - 10 I_{n-1}.$$ 
En d\'eduire une proc\'edure pour calculer $I_n$ pour $n$ quelconque.
\item Utiliser la proc\'edure de~(2) pour calculer en
  flottant $I_5$,  $I_{10}$, $I_{20}$ avec Xcas. Les estimations obtenues vous
  paraissent-elles raisonnables (noter par exemple que, pour tout
  $n$, $0 \leq I_n\leq 1$)?
\item Montrer que l'on peut renverser la proc\'edure 
  ci-dessus:
$$
I_{n-1}=\frac{1}{10n}-\frac{I_n}{10}.
$$
A partir d'une estimation grossi\`ere de $I_{30}$, utiliser cette
proc\'edure inverse pour estimer $I_{20}$. Comparer le r\'esultat avec
la valeur de l'int\'egrale fournie par Xcas.
\end{enumerate}

\paragraph*{Exercice 8}
Pour calculer la valeur
num\'erique de la d\'eriv\'ee d'une fonction $f$ suffisamment
r\'eguli\`ere, 
on peut utiliser l'une des formules suivantes~:
$$f_1=\frac{f(x+h)-f(x-h)}{2h}, \quad f_2=\frac{f(x+h)-f(x)}{h}$$
Mais si $h$ est trop petit, alors $x+h=x$ en flottants, et m\^eme si
$x+h\neq h$, on va introduire des erreurs relatives importantes en
faisant la diff\'erence entre deux termes proches si $h$ est trop
petit.
\begin{enumerate}
\item Justifiez les formules $f_1$ et $f_2$ en faisant un
  d\'eveloppement de Taylor de $f$ en $x$.
\item On suppose que les flottants normalis\'es sont repr\'esent\'es
avec une erreur relative de $\varepsilon$ proche de {\tt1e-15}. 
D\'eterminer en fonction de $h$ l'erreur
absolue pour $f(x+h)-f(x-h)$ et $f(x+h)-f(x)$ en supposant
que $f$ et $f'$ sont de taille proche de 1 au voisinage de $x$.
En d\'eduire l'erreur relative.
\item Pour quelle valeur de $h$ obtient-on une erreur relative
  optimale (on tiendra uniquement compte de la premi\`ere 
puissance de $h$ non nulle dans le reste de la premi\`ere question)~? 
On se limitera donc \`a des valeurs de $h$ 
sup\'erieures ou \'egales.
\item Exp\'erimentez avec des fonctions dont on connait l'expression
de la d\'eriv\'ee, par exemple $\sqrt{\exp(\sin(x))}$ en $x=0$.
\item Proposez une formule approchant la d\'eriv\'ee seconde de $f$
en $x$ faisant intervenir $f(x-h),f(x),f(x+h)$ 
et discutez la valeur de $h$ optimale. Exp\'erimentez.
\end{enumerate}

\paragraph*{Probl\`eme 9 (diff\'erences finies en dimension 1)}

Soit $A$ la matrice tridiagonale de taille $N$ donnée par :
\[A=\left(\begin{array}{ccccc}2&-1&&&\\-1&2&-1&&{\bf 0}\\
    &\ddots&\ddots&\ddots&\\{\bf
      0}&&-1&2&-1\\&&&-1&2 \end{array}\right)\]
\begin{enumerate}
\item
Montrer que 
$$u_k=(\sin(\frac{k\pi}{N+1}),...,
\sin(\frac{kj\pi}{N+1}),...,\sin(\frac{kN\pi}{N+1})) $$
est vecteur propre de $A$ et calculer la valeur propre associ\'ee.
En d\'eduire que si $\tilde{C}$ est une matrice diagonale \`a coefficients
positifs, alors $A+\tilde{C}$ a toutes ses valeurs propres strictement
positives et est inversible.
\item On veut r\'esoudre num\'eriquement 
$$-u'{'}(x)+c(x) u(x)=f(x), u(0)=0, u(1)=0, c(x)\geq 0, \quad x \in ]0,1[$$
On discr\'etise $[0,1]$ en $N+1$ intervalles de m\^eme longueur
$h=1/(N+1)$ et on approche $u(nh)$ par $u_n, n \in [1,N]$ (avec
$u(0)=u_0=0$ et $u(1)=u_{N+1}=0$). On pose $f_n=f(nh),
\vec{f}=(f_1,...,f_n)$, $c_n=c(nh)$ et $C$ la matrice de diagonale
$(c_1,...,c_N)$. D\'eterminer le syst\`eme
lin\'eaire que v\'erifie $\vec{u}=(u_1,...,u_N)$ en fonction de
$A/h^2, C, \vec{f}$.
\item Effectuer la r\'esolution num\'erique par une m\'ethode de
  r\'esolution directe. On pourra utiliser \verb|A:=matrix(N,N);|
pour cr\'eer une matrice dense de taille $N$, puis
\verb|A[0..N-1,0..N-1]:=2.0| pour mettre la diagonale \`a 2, etc.
\item Observez le temps de calcul en fonction de
  $N$, par exemple pour $c(x)=x, f(x)=(1+2x-x^2)e^x+x^2-x$\\
 Repr\'esenter graphiquement l'erreur avec la solution exacte
$u(x)=(1-x)(e^x-1)$.
\item M\'ethodes it\'eratives~: on peut utiliser la m\'ethode de Jacobi pour
obtenir une valeur approch\'ee de $\vec{u}=(u_1,...,u_N)$ plus rapidement lorsque $N$
est grand. On \'ecrit $\tilde{A}=D+(\tilde{A}-D)$ o\`u $D$ est la partie diagonale de
$\tilde{A}$, puis le syst\`eme $\tilde{A}\vec{u}=b$ en rempla\c{c}ant $\vec{u}$ par deux termes
successifs d'une suite r\'ecurrente $v_k$ (\`a valeurs dans $\R^N$)~:
$$ Dv_{k+1}+(\tilde{A}-D)v_k=b \ \Rightarrow \ v_{k+1}=v_k+D^{-1}(b-\tilde{A}v_k) $$
Il peut alors \^etre judicieux d'utiliser une matrice creuse pour
stocker $\tilde{A}$, quel est le nombre d'op\'erations \`a faire par
it\'eration~?\\
On peut justifier la convergence de $v_k$ lorsque $C=0$ (et la vitesse de convergence) 
en regardant la norme
euclidienne de la matrice $I-D^{-1}\tilde{A}$ en se pla\c{c}ant dans la base
orthonorm\'ee de vecteurs propres de $\tilde{A}$. Dans le cas
g\'en\'eral, on peut appliquer le th\'eor\`eme de convergence de la
section sur les m\'ethodes alternatives au pivot du cours. Voir aussi
la m\'ethode de Gauss-Seidel et relaxation.
\end{enumerate}

\end{document}