%% Last modified: <Analyse_calcul_diff16.tex modifié par  Guillaume CONNAN le dimanche 3 avril 2016 à 22h 38min 16s>


\chapterimage{differentiel}
\chaptertext{La notion d'\textit{infiniment petit} a  crée des guerres entre les
  savants pendant deux siècles. Elle n'a été formalisée rigoureusement qu'à la fin du XIX\eme
  siècle.  Cependant,  nous allons  rester  fidèle  au cheminement  de  l'esprit
  humain  et  nous  attacher  à   cette  notion  qui  est  pédagogiquement  plus
  intéressante et se rapproche de l'arithmétique des flottants sur machine Cette
  approche historique  et intuitive du calcul  a été magistralement mise  en \oe
  uvre dans \cite{opac-b1096819}.}

\chapter{Un soupçon de calcul différentiel}


\section{Le calcul infinitésimal}


\subsection{Infiniment petits}

La notion de dérivée est née bien avant celle de limite: dès le XVII\eme \textsc{Fermat}, \textsc{Newton}, \textsc{Leibniz}, Jean \textsc{Bernoulli},
vont développer  le calcul  infinitésimal, avant que \textsc{Bolzano} et \textsc{Weierstra\ss}
ne perfectionnent la notion de limite deux siècles plus tard. 

Mais les approches  des savants du XVII\eme s'adaptent bien  à l'utilisation des
flottants sur machine. 

L'informatique nous incite donc à un retour aux sources historiques.




\begin{remarque}[Analyse non standard et hyperréels]
\ill{Il  est à noter  que dans  les années  1960, un  logicien américain
d'origine allemande (encore un!...), Abraham
\textsc{Robinson},  mit  au  point  l'\textit{analyse non  standard}  en
utilisant les \textit{nombres hyperréels}. 
Il s'agit  de  nos  réels
auxquels on adjoint des nombres \textit{infinitésimaux} qui sont}{robinson}
strictement inférieurs en valeur absolue  à tout nombre réel non nul, et
les   nombres  \textit{infiniment   grands}  dont   les   inverses  sont
infinitésimaux.  Ceci  est fait de  manière très rigoureuse et  permet de
justifier une incroyable intuition que \textsc{Fermat}
avait présenté trois siècles plus tôt.

C'est également assez proche de la définition des flottants selon la norme IEEE 754
\end{remarque}

Pour Jean \textsc{Bernoulli} (1691), \textbf{les infiniment petits sont des quantités
qui  peuvent être  ajoutées à  des quantités  finies sans  changer leurs
valeurs. Les courbes sont des polygones à côtés infiniment courts.} 

Voici  qui s'adapte  à la  notion de  sous-capacité et  au tracé  de courbes  sur
machine.


Voyons  ce que  cela donne  pour les  tangentes à  la  courbe d'équation
$y=x^2$:



%\begin{figure}
\begin{center}
  \includegraphics{limfonc.22}
\end{center}
%\end{figure}


Si $x$ augmente de $\Delta x$, alors $y$ devient $y+\Delta y=(x+\Delta x)^2
=x^2+2x\Delta x+(\Delta x)^2$. Or $y=x^2$.

Pour des valeurs de $\Delta x$ non nulles, cela donne donc:

\[\frac{\Delta y}{\Delta x}= 2x+\Delta x\]


C'est  ici qu'on  a  besoin  de la  notion  intuitive d'infiniment  petits...qui
commence pour nous avec la sous-capacité.


\begin{remarque}
\ill{C'est  Joseph-Louis \textsc{Lagrange}(1736-1813)  ou plutôt  Giuseppe Lodovico
\textsc{Lagrangia},  un  mathématicien  Piémontais d'origine  française,
qui  introduisit  les  termes   et  les  notations  que  nous  utilisons
actuellement:
\begin{quote}
Nous  appellerons  la fonction  $f_x$,  \textit{fonction primitive},  par
rapport aux  fonctions $f'_x$,  $f''_x$, \&c. qui  en dérivent,  et nous
appellerons celles-ci, \textit{fonctions dérivées}, par rapport à celle-là.
\end{quote}
Il a aussi  réfuté la notion d'infiniment  petit et a tenté  de fonder l'analyse
sur les séries de \textsc{Taylor}: nous allons en reparler.
}{lagrange}


\end{remarque}



\subsection{Règles de dérivation}\label{sec::regDeriv}


\subsubsection{Dérivée d'une combinaison linéaire}

Soit  $y(x)=a\cdot u(x)+b\cdot  v(x)$  avec $a$  et  $b$ des  constantes
réelles.
En   posant   comme   précédemment  $y+\Delta(y)=y(x+\Delta   x)$   puis
$u+\Delta(u)=u(x+\Delta   x)$  et   $v+\Delta(v)=v(x+\Delta   x)$,  nous
obtenons $\Delta y=a\cdot \Delta u+b\cdot \Delta v$.
Alors:

\begin{theoreme}[Dérivation d'une combinaison linéaire]
$$y=au+bv     \Longrightarrow     \frac{\mathrm{d}y}{\mathrm{d}x}=a\cdot
\frac{\mathrm{d}u}{\mathrm{d}x}+b\cdot    \frac{\mathrm{d}v}{\mathrm{d}x}
\text{ ou bien } y'=a\cdot u'+b\cdot v'$$
\end{theoreme}


\subsection{Dérivée d'un produit}

\begin{recherche}
  Essayez  de prouver  de  même  une règle  de  dérivation d'un  produit
  $y=u\cdot v$.
\end{recherche}



\begin{theoreme}[Dérivation d'un produit]
$$y=u\cdot v \Longrightarrow \hspace{20em}$$
\end{theoreme}


\subsubsection{Dérivée d'un quotient}


Nous  allons d'abord  établir  un résultat  préliminaire.  Pour $x$  \og
petit\fg{} on \og sent\fg{} que $\frac{1}{1-x}\approx 1$.

Par exemple $\frac{1}{1-10^{-5}}\approx 1,0000100001=1+10^{-5}+10^{-10}$.


Essayons d'être plus précis sur l'approximation. Posons
$$\frac{1}{1-x}=1+\delta$$

et essayons  de déterminer  ce $\delta$ sachant  qu'on peut  le négliger
devant $x$.

Multiplions chaque membre de l'équation par $1-x$. On obtient:


$$1=1-x+\delta(1-x)$$

c'est-à-dire, après simplifications:

$$\delta=x+\delta x$$


Ainsi, comme $\delta x$ est négligeable devant $x$, on peut considérer que
$\delta\approx x$ donc:

$$\frac{1}{1-x}\approx 1+x$$

\begin{remarque}[Série géométrique]
\ill{En 1593, François \textsc{Viète} établit avec une méthode similaire
que, pour tout réel $x$ tel que
$|x|<1$, on a:
$$\frac{1}{1-x}=1+x+x^2+x^3+x^4+x^5+\cdots$$}{Viete}
\end{remarque}



\begin{recherche}
  Vous  pouvez  maintenant  essayer  d'établir une  formule  donnant  la
  dérivée d'un quotient $y=\frac{u}{v}$.

\end{recherche}


\begin{theoreme}[Dérivation d'un quotient]
$$y=\frac{u}{v} \Longrightarrow \hspace{20em}$$
\end{theoreme}



\subsubsection{Fonctions composées}


Soit trois fonctions $f$, $g$ et $h$ telles que $h=f\circ g$.
Ainsi, $h(x)=f \bigl( g(x) \bigr)$.

Posons $g(x+\Delta x)=z+\Delta z$ et $h(x+\Delta x)=y+\Delta y$. 

Alors $y+\Delta y=f \bigl( g(x+\Delta x) \bigr)=f(z+\Delta z)$.
Or    $\frac{\Delta   y}{\Delta   x}=\frac{\Delta    y}{\Delta   z}\cdot
\frac{\Delta z}{\Delta x}$. On obtient donc:


\begin{theoreme}[Dérivation de fonctions composées]
$$\frac{{\rm d}   y}{{\rm d}   x}=\frac{{\rm d}    y}{{\rm d}   z}\cdot
\frac{{\rm  d}  z}{{\rm   d}  x}  \text{  ou  }   h'(x)=f'  \bigl(  g(x)
\bigr)\cdot g'(x)$$
\end{theoreme}


Par  exemple, si $h(x)=\left(3x-5\right)^2$,  alors $y=h(x)=f  \bigl( g(x)
\bigr)$ avec $z=g(x)=3x-5$ et  $f(z)=z^2$. On en déduit que
$h'(x)=\frac{{\rm d}    y}{{\rm d}   z}\cdot
\frac{{\rm d} z}{{\rm d} x}=2z\cdot 3=6(3x-5)$.



\section{Comment calculer un logarithme sur machine..au XVII\eme
  siècle?}


\subsection{Calcul d'une approximation d'une racine carrée par la méthode de \textsc{Héron}}\label{sec::heron}



Le mathématicien \textsc{Héron} d'Alexandrie n'avait pas attendu \textsc{Newton} et le
calcul différentiel  pour trouver  une méthode permettant  de déterminer
une approximation de la racine carrée d'un nombre entier positif puisqu'il a
vécu seize siècles avant Sir Isaac.


Si  $x_n$  est une  approximation  strictement  positive  par défaut  de
$\sqrt{a}$, alors $a/x_n$ est  une approximation par excès de $\sqrt{a}$
(pourquoi~?) et vice-versa.

La    moyenne   arithmétique    de   ces    deux    approximations   est
$\frac{1}{2}\left(x_n+\frac{a}{x_n}\right)$  et constitue  une meilleure
approximation que les deux précédentes.

On  peut  montrer c'est  une  approximation  par  excès (en  développant
$\left(x_n-\sqrt{a}\right)^2$ par exemple).

\begin{pythoncode}
def heron(a, x, n) :
    return x if n == 0 else heron(a, (x + a/x) / 2, n - 1)
\end{pythoncode}

Le résultat est assez bluffant:

\begin{pythoncode}
In [5]: heron(2,1,5)
Out[5]: 1.414213562373095

In [6]: 2**0.5
Out[6]: 1.4142135623730951
\end{pythoncode}


Il  y a  cependant des  discussions au  sujet de  l'exactitude des  calculs dans
$\bbv$: nous en reparlerons dans un prochain chapitre quand nous saurons diviser
deux flottants :-)

\subsection{Calcul d'une approximation de log 2 par la méthode de \textsc{Briggs}}
\label{sec::briggs}

\setindent{-0.75cm}
\begin{flushleft}
\begin{minipage}[t]{3cm}
%\begin{figure}
  \includegraphics[width=3.65cm]{stifel}
%\caption{\foreignlanguage{russian}{Николай Петрович Брусенцов}(1925 -)}
%\end{figure}
\end{minipage}
\end{flushleft}
\setindent{\indentationcours}
\vspace{-3.2cm}
\begin{minipage}[t]{\linewidth}


\end{minipage}


Michael \textsc{Stifel} s'est intéressé dès 1544 aux fonctions vérifiant
$\ell(x\cdot  y)=\ell(x)+\ell(y)$  afin de  simplifier  les calculs  des
astronomes et des navigateurs.  Cependant, cela nécessite de disposer de
tables  \textit{logarithmiques} ($\lambda  o \gamma  o  \varsigma$: mot,
relation,  $\alpha  \rho \iota  \theta  \mu  o  \varsigma$: nombre;  les
logarithmes  sont donc  des  relations entre  les  nombres avant  d'être
l'anagramme d'algorithmes...).

Nous allons nous occuper comme Henry \textsc{Briggs} du logarithme de base
10.

Sachant que $\log_{10}(10)=1$ et  que $\log_{10} \left(a^p\right) =p\cdot \log_{10}(a)$
pour        tout         rationnel        $p$,       on en déduit  alors   que     $\log_{10}
\left(\sqrt{10}\right)=\frac{1}{2}$,                                $\log_{10}
\left(\sqrt{\sqrt{10}}\right)=\frac{1}{4}$, etc.

Nous  pouvons  utiliser  l'algorithme   de  \textsc{Horner}  pour  calculer  les
approximations des racines successives de 10 et on obtient le tableau suivant:



\begin{center}
             \begin{tabular}{lccccc}
	     \hline
             Nombres & 10,0000 & 3,1623 & 1,7783 & 1,3335 & 1,0000\\
             Logarithmes & 1 & 0,5 & 0,25 & 0,125 & 0\\
             \hline
             \end{tabular}
             \end{center}
Le problème, c'est que connaître  le logarithme de 3,1623 n'est pas très
intéressant. On aimerait plutôt connaître ceux des entiers de 1 à 10.

C'est ici  qu'intervient l'anglais  Henry \textsc{Briggs} qui  publia en
1617   ses   premières   tables   de  logarithmes   décimaux   contenant
$\numprint{1000}$ valeurs avec quatorze décimales.

Voyons comment il a procédé pour le calcul de $\log(2)$.

Il a commencé par obtenir des valeurs approchées des nombres suivant:

\begin{center}
\begin{tabular}{cc}
  $\sqrt{10}=10^{\frac{1}{2}}$ & $\sqrt{2}=2^{\frac{1}{2}}$\\
$\sqrt{\sqrt{10}}=10^{\frac{1}{2^2}}$ & $\sqrt{\sqrt{2}}=2^{\frac{1}{2^2}}$\\
$\sqrt{\sqrt{\sqrt{10}}}=10^{\frac{1}{2^3}}$ & $\sqrt{\sqrt{\sqrt{2}}}=2^{\frac{1}{2^3}}$\\
$\vdots$ & $\vdots$\\
$\sqrt{\sqrt{\cdots         \sqrt{10}}}=10^{\frac{1}{2^{54}}}$         &
$\sqrt{\sqrt{\cdots \sqrt{2}}}=2^{\frac{1}{2^{54}}}$\\
\end{tabular}
\end{center}

Il a finalement abouti à:


\[
10^{\frac{1}{2^{54}}}\approx
\underbrace{\numprint{1,00000000000000012781914932003235}}_{1+a} \]

\[ 2^{\frac{1}{2^{54}}}\approx \underbrace{\numprint{1,00000000000000003847739796558310}}_{1+b}
\]


Or  on cherche à  déterminer $x=\log(2)$  mais $x=\log(2)\Leftrightarrow
10^x=2$.

D'après les calculs de \textsc{Briggs}:

\[
1+b=2^{\frac{1}{2^{54}}}=\left(10^x\right)^{\frac{1}{2^{54}}}=(1+a)^x
\]

Un  résultat  classique  que  nous   allons  très  bientôt  étudier  montre  que
$(1+a)^x\approx 1+ax$ pour $a$ assez 
proche de zéro et $x>0$. On en déduit que
$1+b\approx 1+ax$ et donc que:

\[
\log(2)=x\approx
\frac{b}{a}=\frac{\numprint{3847739796558310}}{\numprint{12781914932003235}}
\]

\begin{pythoncode}
In [7]: 3847739796558310 / 12781914932003235
Out[7]: 0.30102999566398114
\end{pythoncode}

Ainsi, $\log(2)\approx \numprint{0,30102999566398114}$.


Cela  nous permet  alors d'avoir  également  $\log(4)=2\times\log(2)$ et
$\log(8)=3\times\log(2)$ mais aussi $\log(5)=\log(10)-\log(2)$.

Il  nous faudrait  donc  obtenir  de même  $\log(3)$  et $\log(7)$  pour
compléter    notre   collection    car    $\log(6)=\log(3)+\log(2)$   et
$\log(9)=2\cdot\log(3)$.

On  peut donc  à titre  d'exercice  calculer ces  deux logarithmes  pour
obtenir:

\[
\log(3)\approx \numprint{0,47712125471966244} \quad \log(7)\approx \numprint{0.8450980400142567}
\]


La   route   est   longue   jusqu'aux   $\numprint{1000}$   logarithmes   de
\textsc{Briggs}\citep{roegel:inria-00543939}:


  \begin{center}
    \includegraphics[height=6cm]{briggs}
  \end{center}

On  peut  retrouver  ces  résultats  sur Sage  qui  permet  de  travailler  très
simplement avec MPFR:

\sagec

\begin{pythoncode}
# Travail commode avec la bibliothèque MPFR de C  via sage
# On ouvre dans emacs un fichier d'extension sage qui met emacs en mode sage
# La syntaxe de programmation est celle de Python
#
 
# On travaille sur 128 bits (on donne la précision = taille de la mantisse)
RN = RealField(113)
 
# Algorithme de Héron pour calculer ?..
def racine(a,x,n):
    if n == 0:
        return RN(x)
    else:
        return racine(RN(a), ??? * RN(0.5), n - 1)
 
# quid ?
def racine_n(a,n):
    if n == 0:
        return RN(a)
    else:
        return racine_n(racine(RN(a),1.0,10), n - 1)
 
# Calcul de log 2
a = ???
b = ???
 
logDeux = RealField(53)(b / a) # résultat projeté sur 64 bits
diffLog = logDeux - (log(2.))/ log(10.) # On compare avec le log 2 de sage
 
# Calcul de log 3 et log 7 ? Généralisation ? log 0 ?  log (-10) ?
\end{pythoncode}

Alors:

\begin{pythoncode}
sage: logDeux
0.301029995663981182
sage: diffLog
5.55111512312578e-17
\end{pythoncode}

\haskell


\subsection{Polynômes interpolateurs}


\subsubsection{Qu'est-ce qu'un polynôme?}

Les  polynômes en  informatique  sont  un outil  très  puissant  et très  utile,
notamment en  théorie de l'information  (détection et correction  d'erreur). Les
entiers en précision  infinie sont considérés comme des  polynômes. Les exemples
abondent mais il s'agit d'une vision discrète  de la notion: un polynôme est une
liste  finie  de  nombres  et  on munit  l'ensemble  des  polynômes  d'opérations
spécifiques.



Nous allons ici étudier plutôt les fonctions polynomiales, i.e. les fonctions de
la forme:

$$p~:~x\mapsto p(x)=a_nx^n+a_{n-1}x^{n-1}+\cdots + a_0$$

où les $a_i$ sont des constantes arbitraires  qui seront pour nous des nombres à
virgule flottante. Si $a_n$ est non nul, on dira que la fonction polynôme est de
degré $n$.

Cette fonction correspond au polynôme:
 $$P=[a_n,a_{n-1},...,a_0]$$


Dans  $\bbf_2$  nous verrons  que  la  distinction  entre polynôme  et  fonction
polynomiale est très  importante. Dans $\bbr$ ou $\bbv$, c'est  moins grave mais
bon...


L'intérêt? Il  s'agit de  fonctions très rapidement  évaluables sur  machine car
faisant  intervenir   uniquement  des  additions   et  des
multiplications.


\subsubsection{Problème d'interpolation}

On  connaît   quelques  valeurs   du  logarithme  donnée   par  la   méthode  de
\textsc{Briggs}.  Cependant, cette  méthode  est très  fastidieuse: il  faudrait
trouver une autre idée pour calculer les valeurs intermédiaires.

Par exemple, on connaît pour l'instant  des approximations des logarithmes de 1,
2, 3 et 4. Comment faire pour déterminer des approximations de tous les nombres entre 1
et 4 avec un pas de 0,25?

Nous allons chercher un polynôme de degré  3 qui passe par ces quatre points. On
cherche donc $a$, $b$, $c$ et $d$ tels que:

$$
\left\{ \begin{array}{l}
p(x)=a+bx+cx^2+dx^3\\
\mathtt{map}\ p\ [1,2,3,4] = \mathtt{map}\ \log\ [1,2,3,4]
\end{array}
\right.
$$


C'est un simple  système linéaire de 4 équations à  quatre inconnues. Cependant,
les équations ont  une forme très spéciale. Voyons une  première éthode proposée
par \textsc{Newton} en 1676.



\subsubsection{Interpolation de \textsc{Newton}}\label{sec::interpolNewton}

Notons $[y_1,y_2,y_3,y_4]= \mathtt{map}\ \log\ [1,2,3,4]$.


\begin{figure}[h!]
  \begin{center}
    \includegraphics[height=7cm]{newton_diff}
    \caption{Problème  d'interpolation   dans  \textit{Methodus  Differentialis}
      publié par \textsc{Newton} en 1676}
    \label{fig:newton_diff}
  \end{center}
\end{figure}

Appliquons l'algorithme  de Gau\ss{}  à ce  système même si  le génie  allemand ne
naîtra qu'un siècle plus tard:

\begin{center}
\begin{tabular}{|l|l|}
  \hline
$x=1$ & $p(x)=a+b+c+d=y_1$\\
$x=2$ & $p(x)=a+2b+4c+8d=y_2$\\
$x=3$ & $p(x)=a+3b+9c+27d=y_3$\\
$x=4$ & $p(x)=a+4b+16c+64d=y_4$\\
\hline
\end{tabular}
\end{center}

En soustrayant les lignes deux par deux, on fait disparaître $a$:

\begin{center}
\begin{tabular}{|l|}
  \hline
 $y_2-y_1=b+3c+7d=\Delta y_1$\\
$y_3-y_2=b+5c+19d=\Delta y_2$\\
$y_4-y_3=b+7c+37d=\Delta y_3$\\
\hline
\end{tabular}
\end{center}

En soustrayant les lignes deux par deux, on fait disparaître $b$:

\begin{center}
\begin{tabular}{|l|}
  \hline
 $\Delta y_2-\Delta y_1=2c+12d=\Delta^2 y_1$\\
 $\Delta y_3-\Delta y_2=2c+18d=\Delta^2 y_2$\\
\hline
\end{tabular}
\end{center}

En soustrayant les lignes deux par deux, on fait disparaître $c$:


$\Delta^2 y_2-\Delta^2 y_1=6d=\Delta^3 y_1$

ce qui donne:

  $d=\frac{1}{6}\Delta^3 y_1$


  $c=\frac{1}{2}\Delta^2
y_1-\Delta^3y_1$


 $b=\Delta y_1-\frac{3}{2}\Delta^2
y_1+3\Delta^3y_1-\frac{7}{6}\Delta^3 y_1=\Delta y_1-\frac{3}{2}\Delta^2y_1+\frac{11}{6}\Delta^3 y_1$


$a=y_1-\Delta y_1+\frac{3}{2}\Delta^2y_1-\frac{11}{6}\Delta^3 y_1-\frac{1}{2}\Delta^2
y_1+\Delta^3y_1-\frac{1}{6}\Delta^3 y_1=y_1-\Delta y_1+\Delta^2 y_1-\Delta^3 y_1$



Finalement:

\[
P(x)=y_1+(x-1)\Delta y_1+\frac{1}{2}(x-1)(x-2)\Delta^2 y_1+\frac{1}{6}(x-1)(x-2)(x-3)\Delta^3 y_1
\]


\textsc{Newton} avait l'habitude de  réunir les résultats dans un schéma
de cette forme:

\begin{center}
\begin{tabular}{cccc}
$\mathbf{y_1}$ &                      &                       &                         \\
               &$\mathbf{\Delta y_1}$ &                       &                         \\
$y_2$          &                      &$\mathbf{\Delta^2 y_1}$&                         \\
               &$\Delta y_2$          &                       & $\mathbf{\Delta^3 y_1}$ \\
$y_3$          &                      &$\Delta^2 y_2$         &                         \\
               &$\Delta y_3$          &                       &                         \\
$y_4$          &                      &                       &                         \\
\end{tabular}
\end{center}


avec $\Delta y_i=y_{i+1}-y_i$, $\Delta^2y_i=\Delta y_{i+1}-\Delta y_i$ et $\Delta^3y_i=\Delta^2 y_{i+1}-\Delta^2 y_i$.




On peut généraliser à des abscisses successives de 1 à $n$:

\begin{theoreme}[Interpolation  de   \textsc{Newton}  passant  par   les  points
  $(1,y_1)$, $(2,y_2)$,..., $(n,y_n)$]\label{th::newton}
\[
p(x)=y_1+(x-1)\Delta y_1+\frac{1}{2!}(x-1)(x-2)\Delta^2 y_1+\cdots+\frac{1}{(n-1)!}(x-1)\cdots(x-(n-1))\Delta^{n-1} y_1
\]
\end{theoreme}

voire à d'autres valeurs d'entiers successifs:

\begin{theoreme}[Interpolation  de   \textsc{Newton}  passant  par   les  points
  $(α,y_1)$, $(α+1,y_2)$,..., $(α+n-1,y_n)$]\label{th::newton_gen}
\[
p(x)=y_1+(x-α)\Delta y_1+\frac{1}{2!}(x-α)(x-α-1)\Delta^2 y_1+\cdots+\frac{1}{(n-1)!}(x-α)\cdots(x-α-(n-2))\Delta^{n-1} y_1
\]
\end{theoreme}


On peut généraliser à des suites quelconques d'abscisses mais c'est un peu plus compliqué.



\section{Polynômes et erreurs}\label{sec::horner}


Nous allons étudier à la \vref{ex::horner} un moyen d'évaluer une fonction
polynomiale en un nombre particulier.

Cependant, nous allons  être à nouveau confrontés au problème  des erreurs suite
au calcul sur les nombres à virgule flottante.

Considérons par exemple la fonction $x\mapsto (x-2)^9$ au voisinage de 2.

Une petite illustration vaut mieux qu'un long discours...Avec le code suivant:

\begin{pythoncode}
comparePoly = 
  plot X11  [ 
    Data2D [Title "Schéma Horner", Color Blue, Style Lines] [] [(x, exHorn x) | x <-  [1.94, 1.940125 .. 2.06]],
    Function2D [Title "Forme factorisée", Color Green] [Range 1.94 2.06, Step 0.000125] (\x -> exFac x)
    ]
\end{pythoncode}


On obtient:


%\begin{figure}
\begin{center}
  \includegraphics[width=0.7\linewidth]{plot_horner_py.png}
\end{center}
%\end{figure}


Bon, évidemment, on est allé chercher  l'exemple qui fâche. Ailleurs, les choses
marchent bien mieux.

Il existe un algorithme  de Horner compensé comme il en existe  un pour la somme
mais c'est un peu  plus compliqué. Les curieux pourront se  référer à un article
de Philippe \textsc{Langlois} et Nicolas \textsc{Louvet} \citep{lalo08}.

\begin{NB}
\illt{Comme nous  ne sommes  qu'en premier cycle,  nous mettrons de  côté à  partir de
maintenant ces problèmes d'erreur dans les évaluations.

Il faudra tout  de même retenir que ces problèmes  existent et sont omniprésents
dans la  machine: est-ce que  mon joueur  est sous le  feu de l'ennemi?  Je dois
savoir si,  en 3D, il  est à gauche  ou à droite d'une  zone ce qui  revient sur
machine  à calculer  un déterminant  dépendant de  3 paramètres,  bref évaluer  un
polynôme de  trois variables. Si  je me trompe dans  mes calculs, quand  je suis
prêt de la zone frontière, mon joueur peut  mourir dans d'atroces souffrances
(sur mon écran seulement, bien sûr).}{call_duty}{0.6}{0.4}

\end{NB}


\section{Série de \textsc{Taylor}}

\subsection{Dérivées d'ordre supérieur}


Nous avons eu  une approche intuitive de la dérivation  jusqu'à maintenant. Elle
nous suffira la plupart du temps. Cependant  il ne faut pas oublier qu'un manque
de rigueur dans l'étude de ces notions peut entraîner de graves erreurs.

Nous rappelons donc  quelques définitions qui nous servirons  à mieux comprendre
certains problèmes.

\begin{definition}[Dérivabilité]
Soit $f$ une fonction numérique définie sur un intervalle I de $\bbr$.

Soit $a\in I$. On dit que $f$  est dérivable en $a$ si, et seulement si,
le rapport:

$$\tau_a(h)=\frac{f(a+h)-f(a)}{h}$$
admet une limite finie lorsque $a$ tend vers 0.

Dans ce cas on note $f'(a)$ cette limite.

\end{definition}


Rappelons que $f'(a)$ est aussi  le coefficient directeur de la tangente
à la courbe représentative de $f$ au point d'abscisse $a$.

\begin{definition}[D\'{e}riv\'{e}e $k^{\grave{e}me}$]

Si $f^{\prime }$ est d\'{e}rivable sur $I$ on note $f^{\prime \prime }$ ou $%
f^{(2)}$ la d\'{e}riv\'{e}e de $f^{\prime }$ et on l'appelle la d\'{e}riv%
\'{e}e seconde de $f.$ De m\^{e}me si $f^{\prime \prime }$ est d\'{e}rivable
sur $I,$ $f^{\prime \prime \prime }$ ou $f^{(3)},$ la d\'{e}riv\'{e}e de $%
f^{\prime \prime },$ est la d\'{e}riv\'{e}e troisi\`{e}me de $f.$ Plus g\'{e}%
n\'{e}ralement on note $f^{(k)}$ la d\'{e}riv\'{e}e $k^{i\grave{e}me}$ de $f$
avec la convention $f^{(0)}=f,$ $f^{(k+1)}=\left( f^{(k)}\right) ^{\prime }.$

\end{definition}


\begin{definition}[Fonctions de classe $\mathcal{C}^{n},$ de classe $\mathcal{C}%
^{\infty }$]

On dit que $f$ est de classe $\mathcal{C}^{n},$ $n\in \mathbb{N},$ sur $I$
ssi $f$ est $n$ fois d\'{e}rivable et ses $n$ d\'{e}riv\'{e}es sont
continues sur $I.$ On remarquera que si $f^{(n+1)}$ existe sur $I$ alors $f$
est de classe $\mathcal{C}^{n}$ sur $I$ qui s'\'{e}crit $f\in \mathcal{C}%
^{n}(I).$ On dit que $f$ est de classe $\mathcal{C}^{\infty }$ sur $I$ pour
exprimer que $f$ est ind\'{e}finiment d\'{e}rivable et, a fortiori, que
toutes ses d\'{e}riv\'{e}es sont continues.

\end{definition}



\subsection{Relations de domination}


\begin{quote}
  {\it Brooks's Law [prov.]

\og Adding manpower to a late software project makes it later\fg{} -- a result
of  the fact that  the expected  advantage from  splitting work  among N
programmers is $O(N)$, but the complexity and
communications cost associated with  coordinating and then merging their
work is $O(N^2)$}
\begin{flushright}
  in \og The New Hacker's Dictionnary\fg{}\\
\href{http://outpost9.com/reference/jargon/jargon\_17.html\#SEC24}{http://outpost9.com/reference/jargon/jargon\_17.html\#SEC24}
\end{flushright}
\end{quote}

Les notations  de \textsc{Landau}(1877-1938) ont  en fait été  crées par
Paul 
\textsc{Bachmann}(1837-1920) en 1894, mais bon, ce sont tous deux des 
mathématiciens allemands.

Par exemple, si l'on considère l'expression: $$f(n)=n+1+\frac{1}{n}+\frac{75}{n^2}-\frac{765}{n^3}+\frac{\cos(12)}{n^{37}}-\frac{\sqrt{\numprint{765481}}}{n^{412}}$$


Quand  $n$  est  \og  grand\fg{}, disons  $\numprint{10000}$,  alors  on
obtient:

$$f(\numprint{10000})=\numprint{10000}+1+0,0001+0,00000000075-0,0000000000000765+\text{
  peanuts}$$

Tous les  termes après  $n$ comptent  pour du beurre  quand $n$  est \og
grand\fg{}: on a expérimenté ceci sur machine avec la notion de \textit{sous-capacité}.



 Donnons une définition pour plus de clarté:

\begin{definition}[\og Grand O\fg{}]
Soit $f$ et $g$ deux fonctions  de $\bbn$ dans $\bbr$. On dit que $f$
est un \og grand O\fg{} de $g$ et on note $f=O(g)$ ou $f(n)=O(g(n))$ si,
et seulement si, il existe une constante strictement positive C telle que $$|f(n)|\leqslant
C|g(n)|$$ pour tout $n\in \bbn$.
\end{definition}



Dans  l'exemple précédent,  $\frac{1}{n} \leqslant  \frac{1}{1}\times 1$
pour tout entier $n$ supérieur à 1
donc $\frac{1}{n}=O(1)$.

De   même,  $\frac{75}{n^2}   \leqslant  \frac{75}{1^2}\times   1$  donc
$\frac{75}{n^2}=O(1)$ mais on peut dire mieux:
 $\frac{75}{n^2}\leqslant  \frac{75}{1}\times\frac{1}{n}$  et  ainsi  on
 prouve que $\frac{75}{n^2}=O \left(\frac{1}{n}\right)$.

En fait, un grand  O de $g$ est une fonction qui  est au maximum majorée
par un multiple de $g$.

Point  de vue  algorithmique, cela  nous donne  un renseignement  sur la
complexité au pire.

Pour  le   tri  fusion,  on  obtient  pour   $n>1$:  $$K(n)\leqslant  4nc
\left(\lfloor                        \log_2(n)\rfloor+1\right)\leqslant
4nc\log_2(n)\left(1+\frac{1}{\log_2(n)}\right) \leqslant 8cn\log_2(n)$$

On en déduit que $K(n)=O(n\log_2(n))$: la complexité est au pire en
$n\log_2(n)$,  c'est ce  qui nous  importe. Le  rôle des  constantes est
accessoire  car  chercher trop  de  précisions  serait illusoire:  l'\og
oubli\fg{}  de ces constantes  correspond en  fait au  différences entre
langages, processeurs, etc.

On peut  cependant faire  mieux avec le  tri fusion  car on a  aussi une
minoration.

C'est le moment d'introduire une nouvelle définition:

\begin{definition}[\og Grand Oméga\fg{}]
  Soit $f$ et $g$ deux fonctions  de $\bbr$ dans lui-même. On dit que $f$
est un \og grand Oméga\fg{} de $g$ et on note $f=\Omega(g)$ ou $f(n)=\Omega(g(n))$ si,
et seulement si, il existe une constante strictement positive C telle que $$|f(n)|\geqslant
C|g(n)|$$ pour tout $n\in \bbn^*$.
\end{definition}

\begin{remarque}
Comme   $\Omega$  est   une  lettre   grecque,  on   peut,   par  esprit
d'unification, parler de \og grand omicron\fg{} au lieu de \og grand O\fg{}...
\end{remarque}

\begin{remarque}
   $f=\Omega(g)\Longleftrightarrow g=O(f)$...
\end{remarque}


On    montre    donc    facilement    pour    le    tri    fusion    que
$K(n)=\Omega(n\log_2(n))$     grâce     à    l'inégalité     $2cn\lfloor
\log_2(n)\rfloor \leqslant K(n)$.


Ainsi   on  a   dans  ce   cas  en   même  temps   $f=O(n\log_2(n))$  et
$f=\Omega(n\log_2(n))$:  c'est  encore  plus  précis et  nous  incite  à
introduire une nouvelle définition:

\begin{definition}[\og Grand Théta\fg{}]
$$f=\Theta(g)\Longleftrightarrow
\begin{cases}
  f=O(g)\\ f=\Omega(g)
\end{cases}
$$
\end{definition}

Le coût de l'algorithme se trouve donc coincé entre deux valeurs de même
ordre.  On peut  ainsi  dire que  la  complexité du  tri  fusion est  en
$n\log_2(n)$ ce qui  est \textit{souvent} mieux que le tri  par insertion dont la
complexité \textit{au  pire} est  en $n^2$. Mais  attention! \textit{Au  pire} ne
signifie pas \textit{toujours}:  si la liste est déjà  triée, le tri par
insertion est plus économique.


Voici maintenant une petite table pour illustrer les différentes classes
de complexité rencontrées habituellement:

\begin{center}
             \begin{tabularx}{\linewidth}{Y{3}}
	     \whline
             \backslashbox{coût}{$n$}&100&$\numprint{1000}$&$10^6$&$10^9$\\
             \whline
             $\log_2(n)$&$\approx7$&$\approx10$&$\approx 20$&$\approx 30$\\
             \hline
             $n\log_2(n)$&$\approx 665$&$\approx \numprint{10000}$&$\approx  2\cdot  10^{7}$&$\approx 3\cdot 10^{10}$\\
             \hline
             $n^2$&$10^4$&$10^6$&$10^{12}$&$10^{18}$\\
             \hline
             $n^3$&$10^6$&$10^9$&$10^{18}$&$10^{27}$\\
             \hline
             $2^n$&$\approx 10^{30}$&$>10^{300}$&$>10^{10^5}$&$>10^{10^8}$\\
             \whline
             \end{tabularx}
             \end{center}

Gardez en tête que l'âge de l'Univers est environ de $10^{18}$ secondes...


Il existe  également deux autres \og comparateurs\fg{}  que l'on utilise
peu  en algorithmique  mais qui  peuvent s'avérer  utiles  dans d'autres
domaines.

\begin{definition}[\og Petit o\fg{}]
$f=o(g)$ si, et seulement  si, pour toute constante positive $\varepsilon$,
il existe un entier $n_0$ tel 
que, pour tout $n \geqslant n_0$, $$\left|f(n)\right| \leqslant \varepsilon\left|g(n)\right|$$
\end{definition}


On dit alors souvent que $f$ est \textit{négligeable} devant $g$.


Contrairement au  grand O, la majoration  doit se faire  quelque soit la
constante $\varepsilon$ et non pas seulement pour une constante arbitraire.

\begin{definition}[Fonctions équivalentes]
  $f(n)\sim g(n)\Longleftrightarrow f(n)=g(n)+o \left(g(n)\right)$
\end{definition}

Voici  quelques   propriétés  fort   utiles  que  nous   démontrerons  à
l'occasion:

\begin{proprietes}[Manipulation des O]
  \begin{itemize}
  \item $O(f)+O(g)=O(f+g)$;
\item $f=O(f)$;
\item $k\cdot O(f)=O(f)$ si $k$ est une constante;
\item $O \left(O(f)\right)=O(f)$;
\item $O(f)\cdot O(g)=O(f\cdot g)$;
\item $O(f\cdot g)=f\cdot O(g)$.
  \end{itemize}
\end{proprietes}

Par exemple, $2n^3+3n^2+5n=O(n^3)+O(n^3)+O(n^3)=O(n^3)$.




\setindent{-0.75cm}
\begin{flushleft}
\begin{minipage}[t]{3cm}
\begin{figure}
  \includegraphics[width=3.65cm]{taylor}
  \caption{Brook \textsc{Taylor}\\ (1685-1731)}
\end{figure}
\end{minipage}
\end{flushleft}
\setindent{\indentationcours}
\vspace{-6.8cm}
\begin{minipage}[t]{\linewidth}


\end{minipage}


\subsection{Polynôme de \textsc{Taylor}}

Brook \textsc{Taylor} fut un savant anglais  qui reprit les travaux de ses aînés
\textsc{Newton},     \textsc{Leibniz},     Jacques     \textsc{Bernoulli}     et
\textsc{Gregory}.

Musicien,  philosophe, peintre,  il est  surtout célèbre  pour les  formules qui
portent son  nom. Il les obtint  en étudiant le problème  de \textsc{Kepler} qui
concerne le calcul de certaines anomalies dans la révolution des planètes.


Son  but  est  encore  une  fois  d'obtenir  une  interpolation  d'une  fonction
compliquée à l'aide d'un polynôme simple.

Il obtient le résultat suivant:

 
\begin{theoreme}[Polynôme de \textsc{Taylor}]
Soit I un  intervalle, soit $f$ une  fonction de Ιvers $\bbr$  de classe $\CR^n$
avec  $n\in\bbn$   et  soit  $a$  un   élément  de  I.  Il   existe  alors  un
\textit{unique} polynôme  $T_{n,a}$ de degré au  plus $n$ dont les  dérivées jusqu'à
l'ordre $n$ coïncident en $a$ avec celles de $f$.
De plus:

$$
T_{n,a}(x) = f(a) + (x-a)f'(a) + \frac{(x-a)^2}{2!}
f''(a)+\frac{(x-a)^3}{3!}f^{(3)}(a)+\cdots +\frac{(x-a)^n}{n!}f^{(n)}(a)
$$
  
\end{theoreme}

Nous admettrons l'unicité  mais \textit{vous vérifierez}
tout de même que pour tout entier $0 \leqslant k \leqslant n$, $$T_{n,a}^{(k)}(a)=f^{(k)}(a)$$



Par exemple, la fonction exponentielle admet comme polynôme de \textsc{Taylor} à
l'ordre 2 en 0:
\begin{eqnarray*}
  \label{eq:1}
  T_{n,0}(x)&=&\exp(0)+(x-0)\exp'(0)+\frac{(x-0)^2}{2!}\exp''(0)\\
        &=&1+x+\frac{x^2}{2}
\end{eqnarray*}



Le  problème  principal  de  sa   formule  est  d'estimer  l'erreur  commise  en
considérant $T_{n,a}$ à la place de $f$:

$$
f(x)       =  T_{n,a}(x)   +R_{n,a}(x)
$$


Voilà pourquoi  on retrouve  son nom  associé à trois  formules qui  donnent une
meilleure idée  de ce reste (\textit{quelle  erreur commet-on?}) et qui  ont été
établies par trois de ses collègues.


Nous ne ferons pas les preuves qui sont assez techniques.

Nous  admettrons   également  que   le  polynôme  de   \textsc{Taylor}  approche
\textit{mieux}\footnote{Ce qui signifie  que $\left|f(x)-T_{n,a}(x)\right| \leqslant
  \left|f(x)-P(x)\right|$ pour tout polynôme P de degré au plus $n$.} $f$
au voisinage de $a$ que tous les autres polynômes de degré au plus $n$.
\subsection{Formule de \textsc{Taylor-Lagrange}}

\begin{theoreme}
 Si $f$ est de classe $\mathcal{C}^{n}$ sur $\left[ a,x%
\right] $ et si $f^{(n+1)}$ existe sur $\left] a,x\right[ $ alors il existe
au moins un r\'{e}el $c\in \left] a,x\right[ $ tel que : 
\begin{equation*}
f(x)=T_{n,a}(x)+\frac{%
\left( x-a\right) ^{n+1}}{(n+1)!}f^{(n+1)}(c)
\end{equation*}%
On dit que l'on a \'{e}crit la formule de Taylor Lagrange \`{a} l'ordre $n.$
\end{theoreme}

Il est souvent pratique de  réécrire cette formule en posant $x=a+h$. Le
nombre $c$ appartient  donc à l'intervalle $]a,a+h[$. Il  existe donc un
réel $\theta\in]0,1[$ tel que:

$$f(a+h)=\sum\limits_{k=0}^{n}\frac{h^{k}}{k!}f^{(k)}(a)+\frac{h^{n+1}}{(n+1)!}%
f^{(n+1)}(a+\theta h)$$


\subsection{Formule de Taylor-Mac Laurin}

Dans la formule de Taylor-Lagrange on remplace  $a$ par $0.$
On obtient : 

\begin{theoreme}
$$f(x)=\sum\limits_{k=0}^{n}\frac{x^{k}}{k!}f^{(k)}(0)+\frac{x^{n+1}}{(n+1)!}%
f^{(n+1)}(\theta x)$$
\end{theoreme}



\subsection{Formule de Taylor avec reste int\'{e}gral}


Si $f$ est de classe $\mathcal{C}^{n+1}$ sur $I$ contenant $a$:

\begin{theoreme}
$$\forall x\in I,f(x)=T_{n,a}(x)+\int\limits_{a}^{x}\frac{\left( x-t\right) ^{n}}{n!}f^{(n+1)}(t)dt$$
\end{theoreme}


\subsection{Formule de Taylor-Young.}

Si $f$ est de classe $\mathcal{C}^{n}$ sur un intervalle $I$ contenant $a$
alors : 


\begin{theoreme}
$$\forall x\in I,f(x)=T_{n,a}(x)+o\left( (x-a)^{n}\right)$$
\end{theoreme}



\section{D\'{e}veloppements limit\'{e}s}


\subsection{Voisinage}

$a\in \mathbb{R},$ nous appellerons voisinage de $a$ toute partie de $%
\mathbb{R}$ contenant un intervalle ouvert contenant $a.$ L'intervalle $%
\left] a-\alpha ,a+\alpha \right[ ,$ avec $\alpha \in \mathbb{R}^{+\ast },$
est un voisinage de $a,$ nous le noterons $\mathcal{V}_{a}$ et $\mathcal{V}%
_{a}-\left\{ a\right\} $ sera not\'{e} $\mathcal{V}_{a}^{\ast }.$ Un
voisinage \`{a} droite de $a$ est un intervalle du type $\left] a,a+\alpha %
\right[ $ et $\left] a-\alpha ,a\right[ $ est un voisinage \`{a} gauche de $%
a.$ Dans ce qui suit nous ne distinguerons pas ces types de voisinages, en
cons\'{e}quence $\mathcal{V}_{a}^{\ast }$ pourra tout aussi bien aussi d\'{e}%
signer un voisinage \`{a} droite ou \`{a} gauche, lors des applications, le
contexte permettra de les distinguer. Un voisinage de $+\infty $ est un
voisinage du type $\left[ A,+\infty \right[ $ avec $A$ r\'{e}el positif
aussi grand que l'on veut que nous pourrons noter $\mathcal{V}_{+\infty
}^{\ast }$. Notez que $\left] -\infty ;-A\right] $ est un voisinage de $-\infty .$


\subsection{Définition}


\begin{definition}
   On dit que $f,$ d\'{e}finie sur $\mathcal{V}_{a},$ admet un
d\'{e}veloppement limit\'{e} d'ordre $n\in \mathbb{N}$ au voisinage de $a$
(on \'{e}crit en abr\'{e}g\'{e}: \og $f$ admet un $DL_{n}V(a)$\fg{}) s'il existe
un polyn\^{o}me $P_{n}\in \mathbb{R}_{n}\left[ X\right] $ (ensemble des polyn%
\^{o}mes de degr\'{e} inf\'{e}rieur ou \'{e}gal \`{a} $n)$ et une fonction $%
\varepsilon $ v\'{e}rifiant : 
\begin{equation*}
\left\{ 
\begin{array}{c}
\forall x\in \mathcal{V}_{a},f(x)=P_{n}(x-a)+(x-a)^{n}\varepsilon (x) \\ 
\varepsilon (a)=0 \\ 
\underset{x\rightarrow a}{\lim }\varepsilon (x)=0%
\end{array}%
\right.
\end{equation*}


c'est-à-dire $f(x)=P_n(x-a)+o \Bigl( (x-a)^n \Bigr)$.

\end{definition}
Il faut bien remarquer que:

$$
q > p \limp (x-a)^q \underset{a}{=}o \Bigl( (x-a)^p \Bigr)
$$

Plus la puissance $p$ est grande, plus $(x-a)^p$ est petit.


Nous n'irons  pas très loin  dans l'étude  des développements limités  mais elle
nous permettra  d'avoir des approximations polynomiales  de nombreuses fonctions
ce  qui  nous  permettra  de  comprendre de  nombreuses  méthodes  employées  en
informatique.




\section{Exemples de résolutions numériques d'équas diff}

\subsection{Ce que l'on dirait si l'on faisait des maths à l'IUT}

Soit $U$ un ouvert de $ℝ×ℝ^n$ et $f\colon U\longrightarrow ℝ^n$ une
application continue sur l'ouvert $U$.

Alors pour $(t_0,x_0)∈U$, on appelle \emph{solution du problème de
\textsc{Cauchy}} :
\[
\begin{cases}
x'(t) = f(t, x(t))\\
x(t_0) = x_0 \qquad {\mathcal{C}}
\end{cases}
\]
tout couple $(I,x)$ où $I$ est un intervalle contenant $t_0$ et $x\colon
I\longrightarrow ℝ^n$ une solution sur $I$ de l'équation différentielle
$x'(t)= f(t, x(t))$ telle que $ x(t_0) = x_0$.

L'équation différentielle qui apparaît dans $(\mathcal{C})$ sera appelée équation
différentielle scalaire si $n=1$, et système différentiel sinon.


Si l'application $f\colon U\longrightarrow ℝ^n$ est de \emph{classe}
$\mathcal{C}^1$ sur l'ouvert $U$, alors, d'après le théorème de
\textsc{Cauchy-Lipschitz},
pour $(t_0,x_0)∈U$, il existe une unique solution maximale au problème
de \textsc{Cauchy} $(\mathcal{C})$.



Dans la plupart des cas, on ne sait pas résoudre explicitement le problème de
Cauchy $(\mathcal{C})$ ; d'où la nécessité de mettre au point des méthodes
numériques de résolution approchée d'un tel problème.

%---------------------------
\subsection{Méthode d'\textsc{Euler}}
%---------------------------

La méthode d'\textsc{Euler} est la méthode la plus simple pour résoudre numériquement une
équation  différentielle.  Son  intérêt  pratique   est  limité  par  sa  faible
précision, c'est pourquoi elle est largement utilisée :-)


L'idée d'\textsc{Euler} consiste à utiliser l'approximation :
\[
x'(t) ≈ \frac{x(t+h) - x(t)}{h}
\quad \text{ pour } h \text{ petit }
\]
Autrement dit,
\[
x(t+h) ≈ x(t) + h ⋅ x'(t) = x(t) + h ⋅ f(t, x(t))
\]
Partant du point $(t_0, x_0)$, on suit alors la droite de pente $f(t_0, x_0)$
sur l'intervalle de temps $[t_0, t_0+h]$.


On pose alors :
\[
\begin{cases}{t_1 = t_0 + h}\\
{x_1 = x_0 + h ⋅ f(t_0, x_0)}
\end{cases}
\]

De nouveau, partant du point $(t_1, x_1)$, on suit alors la droite de pente $f(t_1, x_1)$
sur l'intervalle de temps $[t_1, t_1+h]$, et ainsi de suite.

On construit ainsi une suite de points de la manière suivante :
\[
∀k∈⟦0, N⟧, \quad
\begin{cases}{t_{k+1} = t_k + h}\\
{x_{k+1} = x_k + h ⋅ f(t_k, x_k)}
\end{cases}
\]
La ligne brisée joignant les points $\bigl\{(t_k, x_k) | k∈⟦0,N⟧\bigr\}$ donnera une
solution approchée de notre équation différentielle.

Il faut bien comprendre qu'à chaque étape, on repart dans la direction d'une
solution exacte de l'équation différentielle, mais qui n'est pas celle qui est
solution du problème de \textsc{Cauchy} initial. Sur le graphique ci-contre, on trace une
solution approchée pour la condition initiale $x(0)=1$ et les courbes intégrales
de notre équation qui passent par les points $(t_k, x_k)$.


%\begin{figure}
\begin{center}
   \includegraphics[scale=0.65]{ode_exponentielle2}
\end{center}
%\end{figure}


Pour juger de la qualité d'une méthode (ou schéma) numérique de résolution d'équations
différentielles, il faut prendre en compte plusieurs critères.

\begin{itemize}
    \item
        \index{erreur de consistance}
        L'\emph{erreur de consistance} donne l'ordre de grandeur de l'erreur
        effectué à chaque pas. Par exemple, dans le cas de la méthode
        d'\textsc{Euler},
        l'erreur de consistance mesure l'erreur qu'entraîne le fait d'approcher
        le nombre dérivé par un taux d'accroissement.
        La sommation des erreurs de consistance à chaque pas donnera l'ordre de
        grandeur de l'erreur globale (sous des hypothèses de régularité pour la
        fonction $f$).
        À l'aide de la formule de \textsc{Taylor-Lagrange}, on peut montrer que dans le
        cas de la méthode d'\textsc{Euler}, l'erreur de consistance est dominée par $h^2$
        ; on dit alors que cette méthode est d'ordre 1.
        Ainsi, d'un point de vue théorique,  plus le pas est petit, meilleure
        sera l'approximation.

        \index{schéma numérique!consistant}
        Un schéma numérique est dit \emph{consistant} si la somme des erreurs de
        consistance tend vers 0 quand le pas $h$ tend vers 0 ; ce qui est le cas
        de la méthode d'\textsc{Euler} (en gros, car $Nh^2∼h\to0$ où $N$ désigne le
        nombre d'itérations).

    \item
        \index{schéma numérique!stable}
        La \emph{stabilité} contrôle la différence entre deux solutions
        approchées correspondant à deux conditions initiales voisines : un
        schéma est dit \emph{stable} si un petit écart entre les conditions
        initiales $x(t_0)=x_0$ et $\tilde{x}(t_0)=x_0 + ε$ et de petites erreurs
        d'arrondi dans le calcul récurrent des $\tilde{x}_k$ provoquent une erreur
        finale $\bigl|x_k-\tilde{x}_k\bigr|$ contrôlable.
        La méthode d'\textsc{Euler} est une méthode stable.

    \item
        \index{schéma numérique!convergent}
        Un schéma est dit \emph{convergent} lorsque l'erreur globale (qui
        est le maximum des écarts entre la solution exacte et la solution
        approchée) tend vers 0 lorsque le pas tend vers 0. En fait, pour que le
        schéma soit convergent, il suffit que la méthode soit consistante
        \emph{et} stable.

    \item Enfin pour la mise en application du schéma, il faut aussi prendre en
        compte l'influence des erreurs d'arrondi. En effet, afin de minimiser
        l'erreur globale théorique, on pourrait être tenté d'appliquer la
        méthode d'\textsc{Euler} avec un pas très petit, par exemple de l'ordre de
        $10^{-16}$, mais ce faisant, outre que le temps de calcul deviendrait
        irréaliste, les erreurs d'arrondi feraient diverger la solution
        approchée très rapidement, puisque pour des flottants Python de l'ordre de
        $10^{-16}$, les calculs ne sont plus du tout exacts !

        En pratique, il faut prendre $h$ assez petit (pour que la méthode
        converge assez rapidement du point de vue théorique), mais pas trop
        petit non plus (pour que les erreurs d'arrondis ne donnent pas lieu à
        des résultats incohérents, et que les calculs puissent être effectués en
        un temps fini). La question de trouver de manière théorique
        le pas optimal peut s'avérer un problème épineux. Dans un premier temps,
        on peut se contenter de faire des tests pratiques comme dans l'exemple
        suivant.
\end{itemize}



Dans le script suivant, on applique la méthode d'\textsc{Euler} sur l'intervalle
$[0, 1]$ avec la condition initiale $x(0)=1$.

\begin{pythoncode}
from math import exp


print("{:^10} | {:^12} | {:^12} | {:^13}".format('h','x(t)','exp(t)','erreur'))
print('-'*57)
for i in range(1, 8):
    h, t, x  = 10**(-i), 0, 1
    while t < 1:
        t, x = t + h, x * (h + 1)
    print("{0:>10g} | {1:>.10f} | {2:>.10f} | {3:> .10f}"
            .format(h, x, exp(t),  exp(t) - x))
\end{pythoncode}


\begin{pythoncode}
    h      |     x(t)     |    exp(t)    |    erreur
---------------------------------------------------------
       0.1 | 2.8531167061 | 3.0041660239 |  0.1510493178
      0.01 | 2.7048138294 | 2.7182818285 |  0.0134679990
     0.001 | 2.7169239322 | 2.7182818285 |  0.0013578962
    0.0001 | 2.7184177414 | 2.7185536702 |  0.0001359288
     1e-05 | 2.7182954199 | 2.7183090114 |  0.0000135915
     1e-06 | 2.7182804691 | 2.7182818285 |  0.0000013594
     1e-07 | 2.7182819660 | 2.7182820996 |  0.0000001336
     1e-08 | 2.7182817983 | 2.7182818347 |  0.0000000363
     1e-09 | 2.7182820738 | 2.7182818299 | -0.0000002438
\end{pythoncode}

Pour un pas $h\lesssim 10^{-9}$, le temps de calcul devient trop important, de
plus, on voit poindre l'effet des erreurs d'arrondis, puisque la solution
approchée est devenu supérieure à l'exponentielle, ce qui est mathématiquement
faux, en vertu de l'inégalité :
\[
∀n∈ℕ, \quad \left(1+\frac{1}{n}\right)^n ≤ e
\]

En général, il ne suffit pas qu'un schéma numérique soit convergent pour qu'il
donne de bons résultats sur n'importe quelle équation différentielle.
Encore faut-il que le problème soit mathématiquement bien posé (en
particulier, les hypothèses du théorème de \textsc{Cauchy-Lipschitz} doivent être
vérifiées), qu'il soit numériquement bien posé (continuité suffisamment bonne
par rapport aux conditions initiales) et qu'il soit bien conditionné (temps de
calcul raisonnable).



Considérons par exemple le problème de \textsc{Cauchy} :
\[
\begin{cases}
{x(1)=1}\\
{x'(t) = 3\frac{x(t)}{t} - \frac{5}{t^3}}
\end{cases}
\]
dont la solution est la fonction ???%$φ\colon t\mapsto 1/t^2$.



\begin{center}
    \includegraphics[scale=0.9]{mal_pose}
\end{center}

On constate qu'ici la solution approchée s'écarte assez rapidement de la
solution exacte $φ$. En effet, la forme générale des solutions de l'équation
différentielle précédente est donnée par $t\mapsto λx^3 + 1/x^2$.  La condition
$x(1)=1$ impose $λ=0$, mais dès qu'on s'écarte de la solution $φ$ on suit
\og tangentiellement\fg{} des courbes intégrales qui comportent un terme en
$λx^3$ et donc qui diffère notablement de la solution.  Donc le problème est
\emph{mal posé}. Ici le schéma numérique n'est pas en cause.

Soit à résoudre à présent le problème de \textsc{Cauchy} :
\[
\begin{cases}
{x(1)=1}\\
{x'(t) =100\,(\sin(t)-x)}
\end{cases}
\]
dont la solution est la fonction $φ\colon t\mapsto \frac1{10001} (-100\cos t +
10000\sin t + 100\exp(-100⋅t)) ≈\sin t$.
Que se passe-t-il si on résout ce problème avec un pas de l'ordre de $0.02$ ?


\begin{center}
    \includegraphics[scale=1]{mal_conditionne}
\end{center}
Avec un tel pas, la solution approchée oscille en s'éloignant de plus en plus de la
solution exacte.
En effet, le schéma numérique est donné par :
\[
x_{k+1} = x_k + 100\,h\, \left(\sin t_k - x_k\right)
= (1 - 100h) x_k + 100h \sin t_k
\]
Donc une erreur de $ε_k$ sur $x_k$ aura les répercussions suivantes sur le calcul
des termes ultérieurs :
\[
ε_{k+1} = (1-100h)ε_k \implies ε_{k+n} = (1-100h)^n ε_k
\]
Ainsi donc, tant que $1-100h≥-1$, c'est-à-dire tant que $h\lesssim 0.002$, une petite
erreur sur l'un des termes aura des répercussions allant en s'amplifiant.

Bien que le problème soit numériquement bien posé, il est en fait  \emph{mal
conditionné}.




%---------------------------
\subsection{Méthode de \textsc{Runge-Kutta} d'ordre 4}
%---------------------------



Si, comme nous l'avons déjà dit, la méthode d'\textsc{Euler} présente un intérêt
théorique, on préfère en pratique des méthodes d'ordre plus élevé.
Parmi la multitude des schémas numériques (méthodes à un pas comme celle de
\textsc{Taylor}, méthodes à pas
multiples comme celles d'\textsc{Adams-Bashforth}, d'\textsc{Adams-Moulton}, de
prédiction-correction) celle qui présente le meilleur rapport
précision/complexité est certainement celle de \textsc{Runge-Kutta} d'ordre 4.

En voici le schéma pour $∀k∈⟦0,N⟧$,
\[
\left\lbrace{}
\begin{array}{rcl}
t_{k+1} &=& t_k + h \\
x_{k+1} &=& x_k + \frac16 \left(k_1+2k_2+2k_3+k_4\right)
\end{array}\right.
\quad \text{ où }\quad
\left\lbrace{}
\begin{array}{rcl}
    k_1 &=& h⋅f(t_n, x_n) \\
    k_2 &=& h⋅f\left(t_n + \frac{h}2, x_n + \frac{k_1}2\right) \\
    k_3 &=& h⋅f\left(t_n + \frac{h}2, x_n + \frac{k_2}2\right) \\
    k_4 &=& h⋅f\left(t_n + h, x_n + k_3\right)
\end{array}\right.
\]
Les coefficients qui apparaissent dans ces formules mystérieuses sont
judicieusement ajustés pour obtenir une méthode d'ordre 4, sans pour autant
avoir à calculer les dérivées successives de $f$ (comme c'est le cas dans les
méthodes de \textsc{Taylor}), ou à recourir à une formule de récurrence d'ordre au moins
2 pour définir $x_k$ (comme c'est le cas dans les méthodes à pas multiples).




Définissez  une  fonction  adaptée  à  ce schéma.  Vous  devriez  observer  pour
l'exponentielle


$$
\begin{cases}
{x'(t) = x(t)}\\
{x(0) = 1}
\end{cases}
\qquad{\mathcal{C}}
$$

\begin{pythoncode}
    h      |   erreur euler    |    erreur rk4
----------------------------------------------------
       0.1 |  1.5104931784e-01 |  2.5338874514e-06
      0.01 |  1.3467999038e-02 |  2.2464341498e-10
     0.001 |  1.3578962232e-03 |  2.2204460493e-14
    0.0001 |  1.3592881559e-04 | -2.6645352591e-13
     1e-05 |  1.3591551171e-05 | -5.2180482157e-12
     1e-06 |  1.3591611072e-06 |  2.1464163780e-11
\end{pythoncode}



On constate que pour $h=0.001$, la méthode d'\textsc{Euler} ne donne que trois décimales
significatives du nombre $e$, alors que la méthode de \textsc{Runge-Kutta} en donne
quatorze ! On remarque également qu'il est inutile d'appliquer la méthode de
\textsc{Runge-Kutta} avec un pas $h\gtrsim 0.001$, puisque dans ce cas les erreurs
d'arrondi prennent le dessus par rapport au gain théorique de précision.


\subsection{Système}

Soit à représenter les courbes intégrales du système autonome suivant, appelé
oscillateur de van der Pol :
\[
\begin{cases}
\frac{\dd x}{\dd t} = y\\
\frac{\dd y}{\dd t} = \frac y2 -x -y^3
\end{cases}
\]

Les variables $x$,
$k1$,... $k4$ représentent à présent des listes. Or une instruction comme
$x+k1/2$ renverra un message d'erreur. En effet, la division d'une liste par un
entier n'est pas définie ; de plus, l'addition de deux listes effectue la
concaténation des listes, et non l'addition des listes composante par
composante.

Il faut alors travailler non plus avec des listes, mais avec des objets
supportant la vectorisation, comme c'est le cas des tableaux de la librairie de
tierce partie \verb?NumPy? ; en outre, cette option a l'avantage de réduire de
manière significative les temps de calcul.


%\begin{figure}
\begin{center}
 \includegraphics[scale=.65]{syst_autonome}
\end{center}
%\end{figure}


%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
\subsubsection{Un exemple d'équation différentielle linéaire d'ordre 2}
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

On rappelle qu'une équation différentielle linéaire d'ordre $n$ peut toujours
s'écrire comme un système différentiel d'ordre 1. Illustrons ceci dans le cas
$n=2$ :
en notant :
\[
A(t)  = \begin{pmatrix} 0 & 1 \\ -b(t) & -a(t) \end{pmatrix},\quad
X(t)=\begin{pmatrix} x(t)\\ x'(t) \end{pmatrix},\quad
C(t)  = \begin{pmatrix} 0  \\ c(t)  \end{pmatrix}
\]
on a les équivalences

\begin{align*}
    ∀t∈I\quad x''+a(t)x'+b(t)x = c(t)
& \iff ∀t∈I\quad \begin{cases}
    x' =        x'    \\
    x''= -bx - ax' + c 
  \end{cases}
  \\[1ex]
    &\iff ∀t∈I\quad
    \frac{\dd}{\dd t}\begin{pmatrix} x\\ x' \end{pmatrix}
    = \begin{pmatrix} 0 & 1 \\ -b & -a \end{pmatrix}
    ⋅\begin{pmatrix} x\\ x' \end{pmatrix}
    + \begin{pmatrix} 0  \\ c  \end{pmatrix} \\[1ex]
    &\iff ∀t∈I\quad X' = A(t)⋅X +C(t)
\end{align*}

\index{équation!de \textsc{Bessel} }
Considérons, par exemple,  l'équation différentielle du 2\ieme{} ordre dite de
\textsc{Bessel}  :
\[
t^2 x''(t) + t x'(t) + \left(t^2-α^2\right) x(t) = 0
\iff
X'(t) = A(t)⋅X(t)
\enskip \text{ avec }\enskip
A(t) = \begin{pmatrix}
    0 & 1 \\
    \tfrac{α^2-t^2}{t^2} & -\tfrac{1}{t}
\end{pmatrix}
\]


\begin{center}
    \includegraphics[scale=1]{bessel}
\end{center}



\section{Approximation   de   Pi   et   calcul   approché  d'intégrale au petit bonheur}




On s'intéresse  à l'intégrale sur  $[0,1]$ de la  fonction $f~:~t\mapsto
\sqrt{1-t^2}$: quel est son rapport avec $\pi$?

Nous  utiliserons   Python  pour  rester  basique  et   mieux  voir  les
problèmes. Nous pourrons  cependant aller du côté de  Sage pour avoir un
beau graphique ou faire un calcul en précision infinie.

\sagec

\begin{pythoncode}
sage: P = plot(sqrt(1-x^2),xmin=0,xmax=1)
sage: P.show(aspect_ratio=1)
\end{pythoncode}


\pyth

\subsection{Méthode des rectangles}


C'est a  priori la plus  grossière mais la  plus simple à  écrire. Petit
rappel:

$$\displaystyle\int_{a}^{b}f(x)\;\text{d}x=\displaystyle\sum_{k=0}^{n}\frac{b-a}{n}f
\left(a+k \frac{b-a}{n}\right)+E_n$$

avec $E_n$ l'erreur commise.

% \begin{center}
% \psset{unit=5cm,algebraic}
% \begin{pspicture}(-0.5,-0.5)(1,1.1)
% \psaxes{->}(1.1,1.1)
% \psplot[plotpoints=10,linewidth=1.5pt]{0}{1}{x^2}
% \psStep[linecolor=magenta,StepType=upper](0,1){10}{x^2}
% \psStep[linecolor=blue,fillstyle=vlines](0,1){10}{x^2 }
% \end{pspicture}
% \end{center}


On pourrait avoir cette idée pour calculer cette approximation sur Python:
\pythonc

\begin{pythoncode}
def int_rec_droite(f,a,b,N):
    S = 0
    dt = (b-a)/N
    for k in range(N):
        S += f(a + k*dt)*dt
    return S
\end{pythoncode}

Écrivez la fonction jumelle \verb+int_rec_gauche(f,a,b,N)+.


Il \textit{semble} falloir dix millions d'itérations pour obtenir six bonnes décimales de
$\pi$...

\begin{pythoncode}
In [1]: from math import sqrt,pi

In [4]: for k in range(8): print(pi-4*int_rec_gauche(f,0.,1.,10**k))
3.141592653589793
0.2370743273414746
0.021175621810747725
0.0020371866787654014
0.00020117597846525115
2.0037187828503278e-05
2.001175991139803e-06
2.0003720369032862e-07
\end{pythoncode}



\subsubsection{Méthode des trapèzes}\index{méthode!des trapèzes}


C'est  déjà un  peu plus  précis,  à vue  d'\oe il.  Créer une  fonction
\verb+int_trap(f,a,b,N)+. et faites les mêmes tests.


\begin{pythoncode}
In [6]: for k in range(8): print(pi-4*int_trap(f,0.,1.,10**k))
-0.32250896154796127
-0.01081877967185152
-0.00034420431021464637
-1.0891323113604301e-05
-3.444348353198734e-07
-1.0892040602783482e-08
-3.4457015019029313e-10
-1.162758778150419e-11
\end{pythoncode}

On multiplie notre fonction par 4 et donc son erreur avec. Soyons plus
malins...

\begin{pythoncode}
In [7]: for k in range(8): print(0.25*pi-int_trap(f,0.,1.,10**k))
-0.08062724038699032
-0.00270469491796288
-8.605107755366159e-05
-2.7228307784010752e-06
-8.610870882996835e-08
-2.7230101506958704e-09
-8.614253754757328e-11
-2.9068969453760474e-12

\end{pythoncode}



Pour plus de précision, on pourrait  être tenté de prendre un plus petit
arc  de  cercle, disons  entre  $\frac{\pi}{3}$  et $\frac{\pi}{2}$,  en
utilisant encore $f$. En effet,  la fonction n’étant pas dérivable en 1,
on a un peu peur que  cela crée des perturbations.  Et intuitivement, on
se  dit qu’en  limitant l’intervalle  d’intégration, on  devrait limiter
l’ampleur de l’approximation.



Avec la méthode \lstinline+for+:

\begin{pythoncode}
In [9]: for k in range(8): print(pi/12-(int_trap(f,0,1/2,10**k)-sqrt(3)/8))
-0.005817179530668071
-6.0117319471197916e-05
-6.014041918356305e-07
-6.014064968251631e-09
-6.013917142055902e-11
-6.002975894148221e-13
-8.43769498715119e-15
-4.418687638008123e-14
      \end{pythoncode}
      



Cela semble mieux fonctionner. Bizarre tout de même ce dernier résultat....

On a utilisé des segments  de droite horizontaux, des segments de droite
obliques...Et   si   l’on  utilisait   des   segments   de  parabole   :
intuitivement, cela colle plus à la courbe, cela devrait être plus précis.



\subsubsection{Méthode de \textsc{Simpson}}\index{méthode!de \textsc{Simpson}}


On interpole la courbe par un arc de parabole. On veut que cet arc passe
par les points extrêmes de la courbe et le point d’abscisse le milieu de
l’intervalle. Pour cela, on va déterminer les $c_i$ tels que:
\[\int_a^bf(x) \text{d}x=c_0f(a)+c_1f \left(\frac{a+b}{2}\right)+c_2f(b)\]

soit exacte pour $f(x)$ successivement égale à $1$, $x$ et $x^2$.

Posons $h=b-a$ et ramenons-nous au cas $a=0$.
On obtient le système suivant~:


\[
\left\{\begin{array}{l}
c_0+c_1+c_2=h\\
c_1+2c_2=h\\
c_1+4c_2=\frac{4}{3}h
\end{array}\right.
\]

Résolvez-le et expliquez le résultat:

$$\int_a^bf(x)                     \text{d}x=\frac{b-a}{6}\left(f(a)+4f
  \left(\frac{a+b}{2}\right)+f(b)\right)$$


Créez \verb+int_simps(f,a,b,N)+.

On  subdivise l'intervalle  d'intégration et  on utilise  la  méthode de
\textsc{Simpson} sur chaque subdivision avec notre petit arc de cercle :

\begin{pythoncode}
In [17]: for k in range(8): print(pi/12-(int_simps(f,0,1/2,10**k)-sqrt(3)/8))
5.500957927112582e-05
6.664729046423901e-09
6.684097719755755e-13
3.885780586188048e-16
-1.7208456881689926e-15
-4.907185768843192e-14
4.030109579389318e-14
-6.844913524872709e-12
      \end{pythoncode}
      
Étrange cette précision fluctuante....

On s’aperçoit qu’on est un peu comme un spécialiste de la tectonique des
plaques prenant des photos d’une plage tous les mois à la même heure. Il
peut tout à fait se produire que le montage des photos mette en évidence
que le niveau de la mer descend et il pourra dire à la radio que les gaz
à effet de serre n'ont aucune influence sur la fonte des glaces...

\medskip

Finalement,  il  va  falloir  faire  un peu  plus  de  mathématiques  et
d'informatique et ne
pas se contenter
d'observer des résultats au petit bonheur.




Avant d'attaquer les mathématiques, il faut avoir à l'esprit que le processeur
compte  en base  2 et  nous en  base 10  et que  son \og  zéro\fg{} vaut
environ $2,2\times  10^{-16}$: cela vient du  fait que la  mantisse a 53
bits ce  qui correspond  en gros  à 16 chiffres  décimaux. Il  faut bien
noter que ce n'est pas zéro mais une valeur
nommée communément (pas seulement sur \textsf{Python}{} ) \verb+epsilon+.

Comme c'est l'ordinateur qui va
compter, il faudrait plutôt chercher à le ménager et en tenir compte.

\medskip

Il nous faut  donc quand même regarder sous le  capot pour comprendre la
panne.

Le problème, c'est que nous  ne travaillons pas en précision infinie. En
chargeant   le   module  \verb+sys+,   on   a   accès   à  la   commande
\verb+float_info.epsilon+ qui donne la différence entre \verb+1.0+ et le
nombre en notation scientifique suivant le plus proche:

\begin{pythoncode}
In [18]: from sys import*
In [19]: float_info.epsilon
Out[19]: 2.220446049250313e-16
\end{pythoncode}

Eh  oui,  la  largeur   des  pipe-lines  du  microprocesseur  n'est  pas
infinie. Ses \og réels\fg{} 
admettent des successeurs.

Mais attention, entre 0 et 1, les choses sont différentes:


\begin{pythoncode}
In [20]: e = float_info.epsilon
In [21]: e
Out[21]: 2.220446049250313e-16
In [22]: e/2
Out[22]: 1.1102230246251565e-16
In [23]: 1 + e/2
Out[23]: 1.0
In [24]: 1 + e/2 == 1
Out[24]: True
In [25]: 0 + e/2
Out[25]: 1.1102230246251565e-16
In [26]: e * 1e16
Out[26]: 2.220446049250313
In [27]: -1 - e/2
Out[27]: -1.0
In [28]: 1 - e/2
Out[28]: 0.9999999999999999
\end{pythoncode}

Cet \verb+epsilon+ n'est pas zéro, rappelons-le, mais détermine à partir
de quand deux
flottants seront considérés comme égaux par le système.


Utilisons à nouveau le test \lstinline+pseudo_egal_float+ introduit à la
section précédente.


\begin{pythoncode}
def pseudo_egal_float(a, b):
    return abs(a - b) <= (float_info.epsilon * min(abs(a), abs(b)))
\end{pythoncode}



Si  on  n'y  prête  pas  attention,  on peut  arriver  à  des  résultats
surprenants:

\begin{pythoncode}
In [30]: pseudo_egal_float(0.1 + 0.1 + 0.1 , 0.3)
Out[30]: True
In [31]: 0.1 + 0.1 + 0.1 == 0.3
Out[31]: False
In [32]: 3 * 0.1 == 0.3 
Out[32]: False
In [33]: 4 * 0.1 == 0.4
Out[33]: True
In [34]: 10 + float_info.epsilon - 10
Out[34]: 0.0
In [35]: 1 + float_info.epsilon - 1
Out[35]: 2.220446049250313e-16
In [36]: 10 + 10*float_info.epsilon - 10
Out[36]: 1.7763568394002505e-15
In [37]: x = 0.1
In [38]: 3*x - 0.3
Out[38]: 5.551115123125783e-17
In [39]: 4*x - 0.4 
Out[39]: 0.0
\end{pythoncode}


Rappelons également que le processeur est plus à l'aise avec les puissances de 2
%
car  une multiplication  par 2  ou une  de ses  puissances revient  à un
décalage dans son écriture binaire.


Ici, $0,25=\frac{1}{4}$ en base 10 donc $\frac{1}{100}=0,01$ en base 2.


Pour $0,1=\frac{1}{10}$ en base 10, on obtient $\frac{1}{1010}$ en base
2; posez la division!

% \index{division (en base 2)}
% {\small
% \begin{Verbatim}
%  1      | 1010
%         ------
%  10     | 0,00011
%  100    |
%  1000   |
%  10000  |
% - 1010  |
% ------  |
%    1100 |
% -  1010 |
% ------- |
%      10 |
% \end{Verbatim}
% }
On retrouve au bout d'un moment un reste précédent donc $\frac{1}{10}$ en base 10
admet en base 2 un développement périodique: \verb+0,0 0011 0011 0011...+.


On peut  alors savoir comment ces  nombres sont codés  sur un ordinateur
disposant  de 53  bits pour  la mantisse,  faire de  même pour  $0,3$ et
$3\times 0,1$ et découvrir pourquoi \verb+0.3-3*0.1+ est non nul.

Pour un  exposé brillant et complet  sur le sujet,  étudier les nombreux
documents que William  \textsc{Kahan}, fondateur de la norme  IEEE 754, a
mis en ligne sur sa page web: \href{http://www.eecs.berkeley.edu/~wkahan/}.


Reprenons  l'affinement successif  de  notre subdivision  mais avec  des
nombres de subdivisions égaux à des puissances de 2:


\begin{pythoncode}
In [60]: for k in range(13): print(pi/12 - (int_simps(f,0,1/2,2**k) - sqrt(3)/8))
5.500957927112582e-05
3.935118889741851e-06
2.5687273891294993e-07
1.624742479444663e-08
1.0185862153733183e-09
6.371086991308061e-11
3.9825920339353615e-12
2.490230244234226e-13
1.5376588891058418e-14
6.661338147750939e-16
-5.551115123125783e-17
-1.1657341758564144e-15
-7.771561172376096e-16
      \end{pythoncode}
      


Tout  a l'air  de bien  fonctionner  jusqu'à $2^{10}$  mais ensuite,  on
arrive aux  alentours de  \verb+epsilon+ et cela  commence à  se détraquer
informatiquement.


En  affinant grossièrement  à coups  de  puissances de  10, nous  étions
passés à côté du problème.

Comme quoi, agir intuitivement
en mathématiques ou  en informatique (ici: \og plus  on subdivise petit,
meilleure   sera   la   précision\fg{})   peut   entraîner   de   graves
erreurs...


La fonction  est aussi à  prendre en considération,  puisqu'elle demande
de soustraire à 1 un tout petit nombre:

\begin{pythoncode}
In [62]: f(1e-9)
Out[62]: 1.0
\end{pythoncode}



De plus, cela explique aussi les différences entre la \og méthode
\verb+while+\fg{} et la \og méthode \verb+for+\fg{}.

Le  test   \verb|while  t+dt<=b|   peut  s'arrêter  pour   de  mauvaises
raisons. Par exemple:

\begin{pythoncode}
In [63]: 1 + 0.1 + 0.1 + 0.1 - 0.3 <= 1
Out[63]: False
\end{pythoncode}


La boucle peut  ainsi s'arrêter inopinément, alors qu'on  est loin de la
précision demandée.

%%%%

Précédemment, avec  la méthode \og  gros sabots\fg{}, nous  n'avions pas
vu  de  différence  entre la  méthode  des  trapèzes  et la  méthode  de
\textsc{Simpson} ce  qui contredisait  notre intuition. Observons  à pas
plus feutrés ce qui se passe et incluons les rectangles:


\begin{pythoncode}
In [66]: 
print('-'*75)
print('{:22s} | {:22s} | {:22s}'.format('rectangle gauche', 'trapeze', 'simpson'))
print('-'*75)

for k in range(10):
    print("{:1.16e} | {:1.16e} | {:1.16e}".format(\
        pi/12 - (int_rec_gauche(f, 0, 1/2, 2**k)-sqrt(3)/8), \
        pi/12 - (int_trap(f, 0, 1/2, 2**k)-sqrt(3)/8), \
        pi/12 - (int_simps(f, 0, 1/2, 2**k)-sqrt(3)/8)))

---------------------------------------------------------------------------
rectangle gauche       | trapeze                | simpson               
---------------------------------------------------------------------------
4.5293036853039759e-02 | -5.8171795306680707e-03 | 5.5009579271125819e-05
1.9737928661185844e-02 | -1.4896493887858187e-03 | 3.9351188897418510e-06
9.1241396362000127e-03 | -3.7497837725530836e-04 | 2.5687273891294993e-07
4.3745806294723244e-03 | -9.3912877730750743e-05 | 1.6247424794446630e-08
2.1403338758707591e-03 | -2.3488877122057605e-05 | 1.0185862153733183e-09
1.0584224993743230e-03 | -5.8728876539682062e-06 | 6.3710869913080614e-11
5.2627480586020514e-04 | -1.4682637225482686e-06 | 3.9825920339353615e-12
2.6240327106874517e-04 | -3.6706854433798952e-07 | 2.4902302442342261e-13
1.3101810126225910e-04 | -9.1767299481571030e-08 | 1.5376588891058418e-14
6.5463166981361010e-05 | -2.2941835209344674e-08 | 6.6613381477509392e-16
\end{pythoncode}


On distingue mieux cette fois la hiérarchie des méthodes selon leur efficacité.

Nous  n'avons  pas le  choix:  il est  temps  à  présent d'invoquer  les
mathématiques pour  préciser un  peu ces observations,  même si  cela peut
paraître violent de se lancer dans un raisonnement mathématique avec des
calculs, des théorèmes, sans machine.


En fait, ces trois méthodes sont des cas particuliers d'une même méthode:
on utilise une  interpolation polynomiale de la fonction  $f$ de degré 0
pour les  rectangles, de degré  1 pour les  trapèzes et de degré  2 pour
\textsc{Simpson}.


\section{Méthodes de Newton-Cotes et programmation objet}

%----------------------------------------------------------
% les macros suivantes doivent être utilisées en mode mathématique
%----------------------------------------------------------
\def\DeuxLignes#1#2{\left\lbrace\vcenter{\hbox{#1}\hbox{#2}}\right.}            % 02/02/02
\def\DeuxLignesMaths#1#2{\left\lbrace\vcenter{\hbox{$#1$}\hbox{$#2$}}\right.}   % 02/02/02
\def\DeuxLignesMathsBis#1#2{\left.\vcenter{\hbox{$#1$}\hbox{$#2$}}\right\rbrace}   % 02/02/02

\def\TroisLignes#1#2#3{\left\lbrace\vcenter{\hbox{#1}\hbox{#2}\hbox{#3}}\right.}
\def\TroisLignesMaths#1#2#3{\left\lbrace\vcenter{\hbox{$#1$}\hbox{$#2$}\hbox{$#3$}}\right.} % 22/09/03

% ...ou...
\def\Ou#1#2{\left\bracevert\vcenter{\baselineskip 9pt\hbox{$#1$}\hbox{ou}\hbox{$#2$}}\right.}

\newcommand{\norme}[2][]{\mathopen#1\lVert #2 \mathclose#1\rVert}
  % norme d'une fonction
\newcommand{\norm}{\norme}

\newcommand{\normeinfinie}[2][]{\mathopen#1\lVert #2 \mathclose#1\rVert_\infty}
  % norme d'une fonction

\newcommand{\vabs}[1]{\left\arrowvert   #1 \right\arrowvert}    % valeur absolue
\newcommand{\normi}{\normeinfinie}
\newcommand{\intervalle}[4]{\mathopen{#1}#2 \mathclose{},#3 \mathclose{#4}}
\newcommand{\itvf}[2]{\intervalle{[}{#1}{#2}{]}}
\newcommand{\itvof}[2]{\intervalle{]}{#1}{#2}{]}}
\newcommand{\itvfo}[2]{\intervalle{[}{#1}{#2}{[}}
\newcommand{\itvo}[2]{\intervalle{]}{#1}{#2}{[}}

Voici  une  présentation  plus  professionnelle, aussi  bien  au  niveau
mathématique qu'informatique mais il  aurait été abrupt de commencer par
là...



Pour donner une valeur approchée de l'intégrale $I$, une première idée est de remplacer la
fonction $f$ par un polynôme $P$ qui interpole $f$ en plusieurs points.
Cependant, dès qu'on augmente le nombre de pivots, on risque fort
de se retrouver confronté au phénomène de \textsc{Runge}.

On modifie alors l'idée de départ ainsi : on commence par subdiviser régulièrement
l'intervalle d'intégration en $n$ sous-intervalles $\itvf{x_j}{x_{j+1}}$ où
\[
x_0 = a,\quad x_n = b,\quad \text{et} \quad
\forall j\in \llbracket 0, n\rrbracket,\quad x_j = a + j h,\quad
\text{ avec } h = \frac{b-a}{n}
\]
Ensuite, on remplace $f$ par un polynôme $P_j$ qui interpole $f$ sur chacun des
\og petits\fg{} segments $\itvf{x_j}{x_{j+1}}$.
\begin{itemize}
    \item Si $P_j$ interpole $f$ au point $x_j$, alors le
        graphe de $P_j$ est une droite horizontale : c'est la
        méthode des rectangles.
    \item Si $P_j$ interpole $f$ aux points $x_j$ et $x_{j+1}$, alors le
        graphe de $P_j$ est une droite affine : c'est la méthode des trapèzes.
    \item Si $P_j$ interpole $f$ aux points $x_j$, $\tfrac{x_j+x_{j+1}}{2}$ et $x_{j+1}$, alors le
        graphe de $P_j$ est une parabole : c'est la méthode de \textsc{Simpson} $1/3$.
\end{itemize}
Enfin, on somme les intégrales de chaque polynôme $P_j$ sur $\itvf{x_j}{x_{j+1}}$, et
on obtient une valeur approchée de $I$ :
\[
I = \sum_{j=0}^{n-1} \int_{x_j}^{x_{j+1}} f(x) \text{d} x
= \sum_{j=0}^{n-1} \left(\int_{x_j}^{x_{j+1}} P_j(x) \text{d} x + E_j\right)
\approx \sum_{j=0}^{n-1} \int_{x_j}^{x_{j+1}} P_j(x) \text{d} x
\]
\index{méthode!de \textsc{Simpson}}
Comme ces méthodes ont déjà été introduites à la section \ref{sec::pi}, on ne
s'attarde que sur la méthode de \textsc{Simpson} $1/3$ composée.

On se place sur l'intervalle $\itvf{x_j}{x_{j+1}}$ et on note
$\xi _j=\tfrac{x_j+x_{j+1}}{2}$ ; alors d'après les formules des différences
divisées, le polynôme d'interpolation de $f$
aux points $(x_j, \xi _j, x_{j+1})$ est
donné par
\[
P_j(x) = f[x_j]  + f[x_j, \xi _j] (x-x_j) + f[x_j, \xi _j, x_{j+1}] (x-x_j)(x-\xi _j)
\]
On en déduit, après calculs,
\[
\int_{x_j}^{x_{j+1}} P_j(x) \text{d} x
= \frac{h}{6}\left(f(x_j) + 4 f(\xi _j) + f(x_{j+1})\right)
\]
Par sommation on obtient
\begin{align*}
I
&\approx \frac{b-a}{6n} \left(
f(x_0) + 4 f\left(\tfrac{x_0+x_1}{2}\right) + 2f(x_1)
+ 4 f\left(\tfrac{x_1+x_2}{2}\right) +
\cdots
+ 2f(x_{n-1}) +  4 f\left(\tfrac{x_{n-1}+x_n}{2}\right) + f(x_n)
\right)  \\
&\approx \frac{b-a}{6n} \left(f(x_0) + 2\sum_{j=1}^{n-1} f(x_j)
+ 4\sum_{j=0}^{n-1} f\left(\tfrac{x_{j}+x_{j+1}}{2}\right)
+ f(x_n) \right) \\
&\approx \frac{b-a}{6n} \sum_{k=0}^{2n} w_k f(a_k)
\quad \text{ avec }\quad
a_k = a + k\cdot\frac{b-a}{2n}
\text{ et }
w_k = \TroisLignesMaths{%
1 \text{ si }k=0 \text{ ou }2n}{%
4 \text{ si } 0<k<2n \text{ et $k$ impair}}{%
2 \text{ si } 0<k<2n \text{ et $k$ pair}}
\end{align*}
On peut montrer  que l'erreur commise
est majorée par
\[
\vabs{E}    \leqslant     \frac{1}{2880}\cdot    h^4\cdot    \big\Vert
f^{(4)}\big\Vert\cdot(b-a)
\tag{**}
\]
On dit alors que la méthode est d'ordre 4 ; noter que la méthode est
exacte pour tout polynôme de degré au plus 3.


Pour conclure cette section, revenons un instant sur l'inégalité $(**)$ : cette
majoration de l'erreur possède un intérêt théorique indéniable ; néanmoins, elle
ne permet pas de déterminer directement l'ordre du pas $h$ à choisir pour
obtenir une approximation de $I$ avec une précision donnée.  En effet, comment
déterminer une borne de la dérivée quatrième de $f$ ?  Numériquement,  ce
problème est trop complexe. En pratique, on applique la méthode avec deux pas
différents ($h$ et $2h$ pour minimiser les évaluations de $f$) et on utilise la
différence des deux approximations numériques comme estimation de l'erreur du
moins bon résultat.







\section{Le nombre Pi et les arctangentes}

\subsection{Le nombre Pi et Machin}


\textsf{Python}{} possède un  module de  calcul  en multiprécision  avec les  quatre
opérations arithmétiques de base et la racine carrée:
\index[xcom]{getcontext().prec}

\begin{pythoncode}
In [72]: from decimal import*
# on regle la precision a 30 chiffres
In [73]: getcontext().prec = 30
# on rentre les nombres entre apostrophes, comme des chaines
In [74]: deux=Decimal('2')
In [75]: deux.sqrt()
Out[75]: Decimal('1.41421356237309504880168872421')
In [76]: getcontext().prec = 50
In [77]: deux.sqrt()
Out[77]: Decimal('1.4142135623730950488016887242096980785696718753769')
\end{pythoncode}




Mais ça  ne va  pas nous avancer  à grand  chose puisque pour  avoir 15
bonnes décimales  avec la méthode  de \textsc{Simpson}, il nous  a fallu
dix millions d'itérations.



Faisons un petit détour par l'histoire...



Tout   commence   à  peu   près   en   1671,   quand  l'écossais   James
\textsc{Gregory}  découvre   la  formule
suivante:

\[\Arctan(x)=\sum_{k=0}^{+\infty} \frac{(-1)^kx^{2k+1}}{2k+1}\]


L'inévitable \textsc{Leibniz}  en publie une démonstration  en 1682 dans
son \textit{Acta Eruditorum} où il
est beaucoup question  de ces notions désuètes que  sont la géométrie et
la trigonométrie.

Dans la même  \oe uvre, il démontre ce que  nous appelons aujourd'hui le
critère     spécial    des     séries    alternées:

\medskip
{\itshape
Soit      $(u_n)$      une      suite      réelle      alternée.      Si
$\left(\left|u_n\right|\right)$  est  décroissante  et converge  vers  0
alors:
\begin{itemize}
\item $\displaystyle\sum_{}^{}u_n$ converge;
\item $\left|R_n\right|\leqslant \left|u_{n+1}\right|$ où $(R_n)$ est la
  suite des restes associés à $\displaystyle\sum_{}^{}u_n$.
\end{itemize}
}

\medskip





On peut  même envisager  d'évoquer une démonstration  de ce  théorème en
terminale puisqu'on y parle essentiellement de suites adjacentes...


On en déduit donc que:


\[\left|\Arctan(x)-\sum_{k=0}^n \frac{(-1)^kx^{2k+1}}{2k+1}\right|<\frac{x^{2n+3}}{2n+3}\]

En 1706, l'anglais John \textsc{Machin} en
démontra la fameuse formule:

$$4\,\arctan\frac{1}{5}-\arctan\frac{1}{239}=\frac{\pi}{4}$$

ce qui  lui permit d'obtenir 100  bonnes décimales de  $\pi$ sans l'aide
d'aucune machine mais avec la formule proposée par son voisin écossais.



On appelle polynômes de \textsc{Gregory} les polynômes $G_n(X)=\displaystyle\sum_{k=0}^n \frac{(-1)^kX^{2k+1}}{2k+1}$.


Écrivez  une  fonction  \verb+greg(a,N)+  qui  donne  une  valeur  de
$G_N(a)$.

Pour cela, on aura besoin de travailler en multiprécision.

On crée une fonction \verb+pi_machin+: que fait-elle?


\begin{pythoncode}
def pi_machin(d):
    getcontext().prec = d + 2
    rap = 2*log(5)/log(10)
    N = int(d/rap) + 1
    return 4*(4*greg(5,N) - greg(239,N))
\end{pythoncode}

On récupère sur Sage une approximation de $\pi$.

\sagec
\begin{pythoncode}
sage: pi.n(digits=1000)
\end{pythoncode}

On la copie dans \textsf{Python}{} et on compare:
\pythonc
\begin{pythoncode}
PI =
Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078
164062862089986280348253421170679821480865132823066470938446095505822317253594
081284811174502841027019385211055596446229489549303819644288109756659334461284
756482337867831652712019091456485669234603486104543266482133936072602491412737
245870066063155881748815209209628292540917153643678925903600113305305488204665
213841469519415116094330572703657595919530921861173819326117931051185480744623
799627495673518857527248912279381830119491298336733624406566430860213949463952
247371907021798609437027705392171762931767523846748184676694051320005681271452
635608277857713427577896091736371787214684409012249534301465495853710507922796
892589235420199561121290219608640344181598136297747713099605187072113499999983
729780499510597317328160963185950244594553469083026425223082533446850352619311
881710100031378387528865875332083814206171776691473035982534904287554687311595
62863882353787593751957781857780532171226806613001927876611195909216420199')
\end{pythoncode}




En une seconde et demie, on a mille bonnes décimales de $\pi$...

\begin{pythoncode}
In [91]: PI-pi_machin(1000)
Out[91]: Decimal('-1.5E-1000')
\end{pythoncode}


% On   peut   alors   utiliser    Sage  pour   déterminer   une   fonction
% \verb+mini_greg(a,d)+ avec  $0<a<1$ déterminant une valeur  de $n$ telle
% que $ \frac{a^{2n+3}}{2n+3} \leqslant 10^{-d}$: faites-le!




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% EXERCICES
%
%%%%%


\vfill

\eject













\modeexercice

\section{EXERCICES}

%
\begin{exercice}[Formule de dérivation]
En  utilisant les  résultats de  la section  \vref{sec::regDeriv}, démontrez  la
formule de  dérivation de $x\mapsto  x^n$ et  $x\mapsto 1/x^n$ pour  tout entier
positif $n$.
\end{exercice}



%
\begin{exercice}[Méthode de \textsc{Briggs}]
Inspirez-vous  de  la méthode  présentée  à  la section~\vref{sec::briggs}  pour
calculer une approximation de $\log(3)$ et $\log(7)$.
\end{exercice}


%
\begin{exercice}[Interpolation de \textsc{Newton}(1): le gag]
Déterminez le polynôme  de \textsc{Newton} passant par les  points (1,1), (2,8),
(3,27) et (4,64).  
\end{exercice}




%
\begin{exercice}[Interpolation de \textsc{Newton}(2)]
Déterminez le polynôme de \textsc{Newton} passant par les points (1,5), (2,2), (3,5), (4,2) et (5,2).  
\end{exercice}

%
% \begin{exercice}[Interpolation du logarithme]
% Déterminez le polynôme interpolateur de \textsc{Newton} qui passe par les points
% $(1,\log 1)$, $(2,\log 2)$, $(3,\log 3)$ et $(4,\log 4)$.

% En  utilisant l'algorithme  de  Horner, déduisez-en  les  valeurs du  logarithme
% décimal des nombres de 1 à 4 avec un pas de 0,5.
  
% \end{exercice}



\begin{exercice}
Écrivez la formule du théorème \vref{th::newton_gen} sans ... mais avec $\sum$ et $\prod$  
\end{exercice}


%
\begin{exercice}[Calculs des coefficients des polynômes de \textsc{Newton}]
  Reprenez les calculs effectués dans la section~\vref{sec::interpolNewton} mais
  avec $n$ points d'abscisses 0 à $n-1$.
\end{exercice}



%
\begin{exercice}[Polynôme de \textsc{Newton}: le programme]

Nous voudrions programmer  l'interpolation du logarithme décimal  sur [1,4] pour
obtenir un tracé tel que celui-ci:


%\begin{figure}
\begin{center}
  \includegraphics[height=10cm]{plot_log_py}
\end{center}
%\end{figure}


obtenu      en      utilisant      la      bibliothèque      \texttt{MatplotLib}
(\mintinline{python}{import matplotlib.pyplot as plt}):


\begin{pythoncode}
In [25]: compareTrace(np.log10, 16, 1, 4, 0, 0.7)
\end{pythoncode}


On pourra créer des fonctions intermédiaires:
\begin{enumerate}


\item Déterminez une fonction \lstinline+ddNewton ::  [a] -> [a]+ qui
  renvoie la  liste des $Δ^iy_1$  du \og triangle\fg{} de  \textsc{Newton} (ceux
  en gras...):

  \begin{pythoncode}
In [60]: ddNewton([5,2,5,2,2])
Out[60]: [5, -3, 6, -12, 21]
  \end{pythoncode}

On ourra  par exemple partir d'un  tableau carré où chaque  ligne \texttt{i} est
remplie de \texttt{Y[i]} si on note \texttt{Y} le paramètre de la fonction.



\item   La   formule   du    théorème~\vref{th::newton}   laisse   entrevoir   un
  \texttt{zipWith} (ou braguette avec). Nous disposons déjà de la liste des $Δ^i$.  Il nous reste à
  construire la \textit{base} des polynômes de \textsc{Newton} 
puis à zipper et additionner (bref faire un produit scalaire : une
liste est un vecteur...) pour créer une fonction \lstinline+interpol+:

\begin{pythoncode}
In [61]: interpol(1.5,[1,8,27,64],1)
Out[61]: 3.375
\end{pythoncode}

\item Pour le tracé, expliquez:


\begin{pythoncode}
def compareTrace(f, xmin, xmax, ymin, ymax):
    X  = np.linspace(xmin, xmax, 500)
    Y  = f(X)
    nb_subdivs = xmax - xmin + 1
    subdiv = np.linspace(xmin, xmax, nb_subdivs) 
    fY = f(subdiv)
    coeffs = ddNewton(fY)
    YY = ( lambda x: interpol( x, subdiv, fY, coeffs ) )(X)
    plt.scatter(subdiv, fY, label = 'Poles')
    plt.plot(X, YY , label = 'Interpolation', ls=':')
    plt.plot(X, Y , label = 'Fonction originale',ls='--')
    plt.axis([xmin*1.2, xmax*1.2, ymin, ymax])
    plt.legend()
    plt.title('Interpolation et différences divisées')
    plt.show()
\end{pythoncode}
  

\end{enumerate}


  
\end{exercice}

%
\begin{exercice}[Algorithme de Horner]\label{ex::horner}
  


\ill{La  méthode que nous  allons voir  porte le  nom du  britannique William
George  \textsc{Horner} (1786  - 1837)  mais  en fait  elle fut  publiée
presque  10  ans  auparavant  par  un  horloger  londonien,  Theophilius
\textsc{Holdred} et simultanément par l'italien Paolo \textsc{Ruffini} (1765 - 1822)
mais fut déjà utilisée par  \textsc{Newton} 150 ans auparavant et par le
chinois \textsc{Zhu  Shije} cinq siècles  plus tôt (vers 1300)  et avant
lui  par le Persan  \textsc{Sharaf al-Din  al-Muzaffar ibn  Muhammad ibn
al-Muzaffar al-Tusi} vers (1100) et avant lui par le
Chinois \textsc{Liu  Hui} (vers 200)  révisant un des  résultats présent
dans \textit{Les Neuf
Chapitres sur l'art mathématique} publié avant la naissance de JC...


Soit $P$ un  polynôme. Il s'agit de  calculer l'évaluation de $P$  en une valeur
particulière.



Prenons   l'exemple   de   $P(x)=3x^5-2x^4+7x^3+2x^2+5x-3$.  Le   calcul
classique nécessite 5 additions et 15 multiplications (vérifiez-le!).}{Ruffini}



On peut faire pas mal d'économies de calcul en suivant le schéma suivant :

\begin{center}
\begin{arrayl}
P(x)&=&\ds\underbrace{a_nx^n+\cdots+a_2x^2+a_1x}_{\text{on met $x$ en facteur}}+a_0\\
            &=&\left(\ds\underbrace{a_nx^{n-1}+\cdots+a_2x}_{\text{on met $x$ en facteur}}+a_1\right)x+a_0\\
            &=& \cdots\\
            &=& \left(\cdots(((a_nx+a_{n-1})x+a_{n-2})x+a_{n-3})x+\cdots\right)x+a_0\\
\end{arrayl}
\end{center}


Ici   cela   donne   $P(x)=(((((3x)-2)x+7)x+2)x+5)x-3$  c'est-à-dire   5
multiplications et 5 additions.

 En fait il y a au maximum $2\times$degré
de P opérations (voire moins avec les zéros).



Notez de plus qu'une machine disposant de  la FMA effectue ce calcul encore plus
exactement.

\begin{pythoncode}
In [64]: horn([3,-2,7,2,5,-3], 1)
Out[64]: 12

In [65]: horn([3,-2,7,2,5,-3], 2)
Out[65]: 27
\end{pythoncode}



Reproduisez le graphe de la section~\vref{sec::horner}.

\end{exercice}


%


\begin{exercice}[Méthode d'\textsc{Euler}]
  On ne  dispose que de  la dérivée d'une  fonction: peut-on tracer  malgré cela
  l'allure de sa courbe représentative?

On  créera une  fonction  Haskell permettant  par exemple  de  tracer la  courbe
représentative de  la fonction dont la  dérivée est la fonction  $x\mapsto x$ et
qui s'annule pour $x=1$.

Généralisez votre  fonction Python  pour faire  de même  avec une  fonction qui
vérifie une équation différentielle $$y'=f(x,y)$$

On commencera par créer la liste des coordonnées des points de la courbe approchée.

On l'utilisera ensuite dans la fonction:

\begin{pythoncode}
def affiche_euler(f, x0, y0, xf, h) :
    xs, ys = euler(f, x0, y0, xf, h)
    plt.plot(xs, ys)
    plt.show()
\end{pythoncode}

On obtient par exemple avec:

\begin{pythoncode}
In [151]: affiche_euler(lambda x,y :y, 0,1,2,1e-4)
\end{pythoncode}

la courbe suivante:


%\begin{figure}
\begin{center}
  \includegraphics[height=7cm]{GraphEuler1.pdf}
\end{center}
%\end{figure}



Essayez  d'étendre au  cas  des équations  différentielles d'ordre  2  : $y''  =
f(x,y,y')$.

% Par exemple:


% \begin{pythoncode}
% *AnalyseNumerique> euler2 (\ x y y' -> - y) 0 0 1 10 (2**10)
% \end{pythoncode}

% donne:


% %\begin{figure}
% \begin{center}
%   \includegraphics[height=7cm]{eulerExpoDeux.pdf}
% \end{center}
% %\end{figure}


\end{exercice}





\begin{exercice}[Construction de l'exponentielle]

  \paragraph*{Une équation différentielle}



\textit{L'équation différentielle} $f'=k\,f$ se retrouve dans de
nombreux problèmes : désintégration des noyaux des atomes d'un
corps radioactif, datation au carbone 14, évolution d'une
population où la croissance est proportionnelle au nombre
d'habitants, etc. Le problème est de trouver une fonction la
satisfaisant.

Par exemple, certains phénomènes en mécanique conduisent à étudier
l'équation différentielle $f''=-f$. Nous connaissons au moins deux
fonctions la satisfaisant : cosinus et sinus.



\paragraph*{Construction approchée du graphe d'une solution par la
méthode d'Euler}




Soit $f$ une fonction dérivable sur $\bbr$ vérifiant $f(0)=1$ et,
pour tout  $x$, $f'(x)=f(x)$.  Utilisez la méthode  d'Euler pour  déterminer une
représentation graphique de $f$ sur $[0,2]$.



\paragraph*{Analyse: étude des propriétés mathématiques d'une
solution}


Soit $f$ une fonction dérivable sur $\bbr$ vérifiant $f(0)=1$
et, pour tout $x$, $f'(x)=k\,f(x)$, avec $k\neq 0$.

\begin{enumerate}

\item Montrez que $f'(0)=k$.

\item Soit $y$ un réel fixé et $g_y$ la fonction définie par
$g_y(x)=f(x+y)f(-x)$

\begin{enumerate}

\item Montrez que $g_y$ est dérivable sur $\bbr$ et calculez
$g'_y(x)$.

\item Calculez $g_y(0)$ et déduisez-en que pour tous $x$ et $y$
réels,
$f(x+y)\,f(-x)=f(y)\qquad \mg{(1)}$

\end{enumerate}

\item Montrez alors successivement que:

\begin{enumerate}

\item pour tout réel $x$, $f(x)f(-x)=1$

\item $f$ ne s'annule pas sur $\bbr$

\item pour tous réels $x$ et $y$, $f(x+y)=f(x)\,f(y)$.

\end{enumerate}

\end{enumerate}


\paragraph*{Unicité de la fonction solution}

Peut-on trouver une autre  fonction, $φ$,  distincte de $f$,
et vérifiant les mêmes propriétés que $f$, à savoir~: $φ$  est
une fonction dérivable sur $\bbr$,  vérifiant $φ(0)=1$ et, pour
tout $x$, $φ'(x)=k\,φ(x)$, avec $k\neq 0$?

\medskip

Comme $f$ ne s'annule pas sur $\bbr$, on peut définir la
fonction $ψ=φ/f$.

Vérifiez que $ψ$ est dérivable sur $\bbr$, calculez sa dérivée.
Que peut-on en déduire pour $ψ$? Montrez alors que $f=φ$.


\paragraph*{Synthèse}

Nous avons cherché des solutions au problème:

\begin{quote}
{\itshape
$f$ une fonction dérivable sur $\bbr$ vérifiant $f(0)=1$ et,
pour tout $x$, $f'(x)=k\,f(x)$, avec $k\neq 0$.
}

\end{quote}


Nous avons montré que, \textit{si une telle fonction existe} (ce
que nous prouverons dans un prochain chapitre), alors elle est
unique et elle  vérifie nécessairement la relation
$f(x+y)=f(x)\,f(y) \qquad (2)$.

Il reste à vérifier que, réciproquement, une fonction dérivable,
non nulle, vérifiant la relation (2)  est nécessairement telle que
$f(0)=1$ et vérifie pour tout réel $x$ $f'(x)=kf(x)$, avec $k$ un
réel non nul.

\bigskip

Cette vérification n'est pas anodine et conclut notre raisonnement
d'\textit{analyse-synthèse}.
%Par exemple, on cherche un plombier
%pouvant se glisser dans les canalisations d'une maison pour la
%réparer.
%
%Analyse : nécessairement, cet homme devra mesurer 1 cm de haut,
%avoir des branchies pour respirer sous l'eau, être nyctalope.
%Synthèse :


\begin{enumerate}

\item Montrez que $f$ ne s'annule pas et que $f$ est à valeurs
strictement positives.

\item Montrez que, comme $f$ n'est pas la fonction nulle, alors
$f(0)=1$ en utilisant la relation (2).

\item Soit $a$ un réel fixé. On définit la fonction $φ~:~x\mapsto
f(x+a)$ et la fonction $ψ~:~x\mapsto f(x)\times f(a)$.

Montrez que $f'(x+a)=f(a)\times f'(x)$, puis que, pour tout réel
$a$, $f'(a)=k\,f(a)$, où $k$ est un réel que vous déterminerez.

\end{enumerate}


\end{exercice}










\begin{exercice}[Lynx et lapins]


  \begin{multicols}{2}
    
\begin{figure}
\begin{center}
  \includegraphics[height=4cm]{lokta}
\end{center}
\caption{Alfred \textsc{Lokta}\\US (1880-1949)}
\end{figure}

\begin{figure}
\begin{center}
  \includegraphics[height=4cm]{volterra}
\end{center}
\caption{Vito Volterra\\It. (1860-1940)}
\end{figure}

  \end{multicols}

Vous connaissez ce problème par c\oe ur...



$$    \left\{\begin{array}{ccc}\dfrac{\mathrm{d}x(t)}{\mathrm{d}t} &=& x(t)\ \Big(\alpha - \beta y(t)\Big) \\ \dfrac{\mathrm{d}y(t)}{\mathrm{d}t} &=& -y(t) \ \Big(\gamma - \delta x(t)\Big) \end{array}\right. 
$$

\begin{itemize}
\item     $t$ est le temps ;
   \item  $x(t)$ est l'effectif des proies en fonction du temps ;
  \item   $y(t)$ est l'effectif des prédateurs en fonction du temps ;
    \item        les       dérivées        $\mathrm{d}x(t)/\mathrm{d}t$       et
      $\mathrm{d}y(t)/\mathrm{d}t$ représentent la  variation des populations au
      cours du temps.

    \end{itemize}
    

Les paramètres suivants caractérisent les interactions entre les deux espèces :
\begin{itemize}
\item     $\alpha$, taux de reproduction des proies (constant, indépendant du nombre de prédateurs) ;
  \item     $\beta$, taux de mortalité des proies dû aux prédateurs rencontrés ;
  \item   $\gamma$, taux de mortalité des prédateurs (constant, indépendant du nombre de proies) ;
   \item  $\delta$, taux de reproduction des prédateurs en fonction des proies rencontrées et mangées ;
\end{itemize}



% \begin{pythoncode}
% def voltera(t,V):
%     x,y = V[0], V[1]
%     return np.array([x*(3 - 2*y), y*(-4 + x)]) 
% \end{pythoncode}


%Prenez par exemple $α=3, β=-2, γ=-1, δ=4$.

Vous connaissez la méthode d'Euler par c\oe ur...

% \begin{pythoncode}
% def EulerV(F,V0,h,n):
%     t = 0
%     V = V0
%     X , Y , T = [V[0]] , [V[1]], [0]
%     for i in range(n):
%         V = V + h*F(t, V)
%         t += h
%         T.append(t)
%         X.append(V[0])
%         Y.append(V[1])
%     plt.plot(T,X)
%     plt.plot(T,Y)
%     # plt.plot(X,Y)
% \end{pythoncode}

Tracez $x$ et $y$ en fonction de $t$ puis $y$ en fonction de $x$.

\end{exercice}



% \begin{exercice}[Deux corps dans l'espace...]
  

% Bla bla Principe Fondamental de la Dynamique bla bla bla:

% $$
% r''-rθ'^2=-\frac{GMt}{r^2}+\frac{Ar}{M_0-At}\qquad rθ''+2r'θ'=\frac{Arθ'}{M_0-At}
% $$


% Mais  comment faire!?  Mais  comment  appliquer la  méthode  d'Euler?  Y a  deux
% inconnues et en plus y a des dérivées secondes...


% On pose $u_1=r$, $u_2=r'$, $u_3=θ$, $u_4=θ'$ et tout s'arrange.


% \end{exercice}



% \begin{exercice}[Centrale 2015 avec une solution exacte du problème des 3 corps]
%  \textit{Sur une idée de \href{http://rouxph.blogspot.fr/}{Philippe \textsc{Roux}}}

% % Voici le PFD appliqué à trois corps célestes:

% % $${d^2\over dt^2}X_i =\sum_{j\neq i}-{{\cal G}m_j(X_i-X_j)\over |X_i-X_j|^3}=\sum_{j\neq i} F_g(X_i-X_j,m_j) ~~~~i=1,2,3$$
% % avec 6 conditions initiales $X_i(0), {d\over dt}X_i(0)$. 

% % $$
% % U=\begin{pmatrix}X_1\\ {d\over dt}X_1\\ X_2\\ {d\over dt}X_2\\ X_3\\{d\over dt}X_3 \end{pmatrix} 
% % ~~~~{d\over dt}U =F(U)=\begin{pmatrix} 
% % U_2\\ 
% % F_g(U_1-U_3,m_2)+F_g(U_1-U_5,m_3)\\ 
% % U_4\\ 
% % F_g(U_3-U_1,m_1)+F_g(U_3-U_5,m_3)\\ 
% % U_6\\ 
% % F_g(U_5-U_3,m_2)+F_g(U_5-U_1,m_1)\\ 
% %  \end{pmatrix}
% % $$

% % Alors:


% % $${d\over dt}U_{2i-1}= {d\over dt}X_i=U_{2i}\qquad {d\over dt}U_{2i}= {d^2\over dt^2}X_i=F_{G,i}(X_1,\dots,X_N)$$

% % Dans   cet  \href{http://arxiv.org/abs/math/0011268}{article},   des  conditions
% % initiales permettent  d'obtenir une solution  périodique pour des corps  de même
% % masse.

% % Si on se place  dans le cas d'un mouvement plan les termes  $U_i$ sont des vecteurs
% % de $\bbr^2$ contenant les positions $x_i,y_i$ et les vitesses $vx_i,vy_i$ des 3 corps :

% % $$
% % U_1=\begin{pmatrix}x_1\\y_1\end{pmatrix},~~U_3=\begin{pmatrix}x_2\\y_2\end{pmatrix},~~U_5=\begin{pmatrix}x_3\\y_3\end{pmatrix}~~~~  
% % U_2=\begin{pmatrix}vx_1\\vy_1\end{pmatrix},~~U_4=\begin{pmatrix}vx_2\\vy_2\end{pmatrix},~~U_6=\begin{pmatrix}vx_3\\vy_3\end{pmatrix}
% % $$



% Nous allons terminer la résolution du problème des N corps de Centrale 2015 avec
% $N=3$ et les conditions initiales suivantes:

% $$
% x_1=−x_2=0.97000436~~~y_1=-y_2=−0.24308753,~~~x3=y_3=0
% $$

% $$
% vx_3=−2vx_1=−2vx_2=−0.93240737~~~vy_3=−2vy_1=−2vy_2=−0.86473146
% $$


% avec des unités telles que les masses et la constante $\cal G$ valent 1.

% %A-t-on un schéma périodique? Quelle est la période?

% On rappelle que la PFD donne $$\ve{F}_{k/j}=G \frac{m_jm_k}{r_{jk}^3}\ve{P_jP_k}$$

% avec $r_{ij}$ la distance entre les corps $j$ et $k$.

% On  utilisera  des  tableaux  de  \verb+numpy+  pour  effectuer  des  opérations
% vectorielles.

% \begin{pythoncode}
% import numpy as np
% import matplotlib.pyplot as plt
% import matplotlib.animation as animation

% G = 1
% m = 1
% position = np.array([[0.97000436,-0.24308753],[-0.97000436,0.24308753],[0,0]])
% #position2 = np.array([[0.99000436,-0.24308753],[-0.97000436,0.24308753],[0,0]])
% vitesse = np.array([[ 0.46620368,  0.43236573],[ 0.46620368,  0.43236573],[-0.93240737,-0.86473146]])

% def norme(v) :
%     """ norme d'un vecteur"""
%     return 

% def force2(m1,p1, m2, p2) :
%     """ force exercée par le corps 2 sur le corps 1"""
%     p1p2 =  
%     return 

% def force(j,m,pos) :
%     """ somme des forces exercées sur le corps j
%         on pourra utiliser np.sum(array, axis = 0)"""
%     return

% def gravit(pos) :
%     """renvoie le vecteur dont les composantes sont les accélérations de chaque corps pour une position donnée """
%     return 

% def verlet(f, pos0, vit0, h, n) :
%     """ on utilise le schéma de Verlet vu en cours vectorialisé"""

% def phase(f, pos0, vit0, h, n, no) :
%     """ trace les mouvements des trois corps"""
% \end{pythoncode}


% On obtient une courbe fermée:


% %\begin{figure}
% \begin{center}
%   \includegraphics[height=5cm]{TroisCorps0}
% \end{center}
% %\end{figure}

% qui ne l'est plus si une valeur initiale est légèrement perturbée (1\%):


% %\begin{figure}
% \begin{center}
%   \includegraphics[height=5cm]{TroisCorps1}
% \end{center}
% %\end{figure}

% et si une valeur initiale est un peu plus perturbée (10\%):


% %\begin{figure}
% \begin{center}
%   \includegraphics[height=5cm]{TroisCorps2}
% \end{center}
% %\end{figure}


% \end{exercice}





\begin{exercice}[MINES-PONTS 2015]


\section*{Introduction}

Les imprimantes sont des systèmes mécatroniques (combinaison synergique et systémique de la mécanique, de l'électronique et de l'informatique en temps réel) fabriqués en grande série dans des usines robotisées. Pour améliorer la qualité des produits vendus, il a été mis en place différents tests de fin de chaîne pour valider l'assemblage  des produits. Pour un de ces tests, un opérateur connecte l'outil de test sur la commande du moteur de déplacement de la tête d'impression et sur la commande du moteur d'avance papier. Une autre connexion permet de récupérer les signaux issus des capteurs de position. Différentes commandes et mesures sont alors exécutées. Ces mesures sont transmises sous la forme d'une suite de caractères vers un ordinateur. Cet ordinateur va effectuer différentes mesures pour valider le fonctionnement de l'électromécanique de l'imprimante. L'ensemble des mesures et des analyses est sauvegardé dans un fichier texte. Afin de minimiser l'espace occupé, les fichiers sont compressés. Une base de données stocke les informations concernant les mesures et les imprimantes sur lesquelles elles portent, et permet à l'entreprise d'améliorer la qualité de la production après diverses études statistiques.
\medbreak

Le problème comporte trois grandes parties indépendantes, sauf par le thème traité.
\begin{itemize}
\item Dans la partie~\ref{donnees}, on s'intéresse au processus de réception par l'ordinateur des données transmises par le capteur.
\item Dans la partie~\ref{numerique}, on procède au traitement des données numériques récupérées à fin de validation de l'imprimante assemblée.
\item Dans la partie~\ref{SQL}, on étudie quelques aspects de la base de données constituée.
\end{itemize}


\section{Réception des données issues de la carte d'acquisition}\label{donnees}
\subsection{Le capteur}
Le capteur utilisé est analogique, et comporte un convertisseur (appelé convertisseur analogique numérique) qui traduit les mesures effectuées sous forme numérique. Chaque donnée est ainsi codée comme un entier.

\begin{ques}
Rappeler un principe de représentation en machine des entiers signés sur $10$ bits. Préci\-ser la plage des valeurs entières représentables selon ce principe.
\end{ques}

\begin{ques}
On considère que les valeurs analogiques s'étendent en pleine échelle de $-5$~V à $5$~V et sont converties en entiers signés représentés sur $10$ bits. Donner une valeur approchée de la résolution de la mesure en volts.
\end{ques}

Dans la suite du problème, la question du codage des entiers en binaire n'intervient plus.

\subsection{Liaison avec l'ordinateur}
Une liaison série asynchrone permet la communication entre la carte de commande/acquisition et le~PC. Les informations correspondant à une mesure sont envoyées par la carte électronique sous la forme d'une suite de caractères, appelée trame. Plusieurs trames consécutives correspondant à différents jeux de mesures peuvent \^etre transmises sans interruption, le format m\^eme des trames devant permettre de les délimiter. Un exemple de trame est :
\begin{center}
\ttfamily
\small
\def\espace{\hskip1.95pt}
\begin{tabular}{|@{\espace}*{20}{@{\espace}c@{\espace}|}}
\hline
'U'&'0'&'0'&'3'&'+'&'0'&'1'&'2'&'+'&'0'&'0'&'4'&'-'&'0'&'2'&'3'&
'9'&'9'&'9'&'3'\\
\hline
\end{tabular}
\end{center}

Une trame est donc constituée de la manière suivante :
\begin{itemize}
\item un caractère d'en-tête, qui est une des trois lettres \texttt{'U'}, \texttt{'I'}, \texttt{'P'}, permettant d'identifier la nature de la mesure effectuée (\texttt{'U'} : tension moteur,  \texttt{'I'} : courant moteur, \texttt{'P'} : position absolue) ;
\item trois chiffres constituant une valeur entière précisant le nombre $N$ de mesures qui constituent la suite de la trame ;
\item un ensemble de $N$ blocs consécutifs de quatre caractères. Chacun des blocs est constitué d'un caractère de signe \texttt{'+'} ou \texttt{'-'}, suivi de trois chiffres donnant une valeur absolue. L'entier signé ainsi codé sur quatre caractères donne une valeur issue de la conversion analogique/numérique d'une mesure ;
\item un dernier bloc de quatre chiffres constituant ce qu'on appelle une somme de contrôle (\textit{checksum}). Celle-ci est calculée en formant la somme des valeurs des données précédentes de la trame, et en prenant le reste de la division euclidienne par $10000$ de cette somme. 
\end{itemize}

La somme de contrôle (\textit{checksum}) permet de vérifier si la trame a été correctement transmise. On peut, à réception d'une trame, calculer la somme de contrôle associée aux données reçues, et la comparer à la somme de contrôle contenue dans la trame. Si ces deux valeurs ne coïncident pas, les données ont été altérées lors de la transmission, et il convient de ne pas les exploiter.

\medbreak

La fonction \texttt{car\_read(nbre\_carac)} permet d'obtenir des caractères consécutifs issus de la liaison asynchrone. Elle prend en argument un entier \texttt{nbre\_carac}, et renvoie une chaîne de caractères constituée de \texttt{nbre\_carac} caractères. Ainsi, dans l'exemple ci-dessus, deux appels consécutifs à \texttt{car\_read(5)} puis \texttt{car\_read(3)}, renverront respectivement les chaî\-nes \texttt{'U003+'} et \texttt{'012'}.

\begin{ques}
Toujours à partir de l'exemple donné plus haut, on suppose que les deux appels \begin{center}
\texttt{car\_read(5)} puis~\texttt{car\_read(3)}
\end{center} ont été effectués. Quelle valeur est alors renvoyée par l'appel~\texttt{car\_read(4)} ?
\end{ques}

Le programme de lecture des données transmises à l'ordinateur lit des caractères jusqu'à obtention d'un caractère~\texttt{x} d'en-tête égal à~\texttt{'U'}, \texttt{'I'} ou~\texttt{'P'}. Il fait alors appel à une fonction de lecture de trame~\texttt{lecture\_trame(x)} lisant les données numériques contenues dans la trame qui commence par~\texttt{x}.

\begin{ques}
Écrire la fonction~\texttt{lecture\_trame(x)} prenant en argument le caractère d'en-tête~\texttt{x} qui vient d'être lu, et renvoyant une liste~\texttt{[x,L,S]}, où \texttt{L} est une liste d'entiers représentant les mesures transmises dans la trame commençant par~\texttt{x}, et~\texttt{S} est la valeur de la somme de contrôle transmise dans la trame commençant par~\texttt{x}.
\end{ques}

Dans la question suivante, on utilise une valeur issue d'un appel à la fonction \begin{center}
\texttt{lecture\_trame(x)}.
\end{center} Une telle valeur sera stockée dans une variable \texttt{trame} \textit{via} une affectation \begin{center}
\texttt{trame = lecture\_trame(x)}.
\end{center}

\begin{ques}
\'Ecrire une fonction \texttt{checkSum(trame)} prenant en argument une liste \texttt{trame}, et renvoyant un booléen indiquant si le troisième élément de la liste \texttt{trame} est égal à la somme de contrôle calculée à partir des valeurs stockées dans le deuxième élément de cette liste.
\end{ques}

\begin{ques}
\'Ecrire une fonction \texttt{lecture\_intensite()}, sans argument, et procédant au traitement suivant :
\begin{itemize}
\item[\ding{172}] obtenir un caractère jusqu'à ce que ce soit le caractère \texttt{'I'} ;
\item[\ding{173}] lire la trame dont on vient de lire le caractère d'en-tête et stocker les informations dans une variable \texttt{trame} ;
\item[\ding{174}] si la  variable \texttt{trame} ne présente pas une somme de contrôle valide, recommencer à~\ding{172};
\item[\ding{175}] renvoyer la liste des mesures numériques alors contenues dans la variable \texttt{trame}.
\end{itemize}
On fait l'hypothèse que le canal de transmission est de qualité suffisante pour assurer qu'on finira par obtenir une trame dont la somme de contrôle est correcte.
\end{ques}

\section{Analyse des mesures}\label{numerique}

{\bfseries Les fonctions demandées pourront être écrites dans le langage~\texttt{Python} (éventuellement à l'aide de sa bibliothèque~\texttt{numpy}) ou dans le langage~\texttt{Scilab}.

On exclut dans cette partie tout appel à une fonction de calcul approché d'intégrale ou de résolution approchée d'équations différentielles qui serait préprogrammée dans le langage choisi ou dans l'une de ses bibliothèques.}

\subsection{Traitement numérique}
\begin{ques}
Rappeler le principe de la méthode des trapèzes permettant de calculer une valeur approchée de l'intégrale d'une fonction $f$ sur un intervalle $[a,b]$. On ne demande pas ici de coder de fonction informatique permettant ce calcul.
\end{ques}

Pour valider le fonctionnement de l'imprimante, on considère une suite de valeurs de mesures de l'intensité du courant (en ampères) du moteur de la tête d'impression, stockée dans un tableau \texttt{mesures}. On s'intéresse alors à la valeur moyenne et à l'écart-type de cette distribution de valeurs.


Pour une fonction $f$ définie (et continue) sur un segment $[0,t_{\text{final}}]$, la valeur moyenne de la fonction sur ce segment est définie par :
\[f_{\text{moy}} = \frac{1}{t_{\text{final}}} \int_0^{t_{\text{final}}} f(t)\mathrm{d}t.\]

On ne dispose ici des valeurs de l'intensité qu'en des instants prédéfinis : les mesures ont été effectuées toutes les $2$ millisecondes.

\begin{ques}
\'Ecrire une fonction \texttt{trapezes(L)} prenant en argument un tableau~\texttt{L} de valeurs numériques, et renvoyant une valeur approchée, calculée par la méthode des trapèzes, de la valeur moyenne d'une grandeur variant durant un certain intervalle de temps $[0,t_{\text{final}}]$, en supposant que le tableau~\texttt{L} contient les valeurs de cette grandeur mesurées toutes les $2$ millisecondes.
\end{ques}

On définit l'écart-type d'une fonction sur un segment $[0,t_{\text{final}}]$ par :
\[f_{\text{ec}} = \sqrt{\frac{1}{t_{\text{final}}}\int_0^{t_{\text{final}}} (f(t) - f_{\text{moy}})^2 \mathrm{d}t}.\]

\begin{ques}
\'Ecrire une fonction \texttt{indicateurs(mesures)} prenant en argument un tableau \texttt{mesures} de flottants codant les valeurs de l'intensité mesurées toutes les $2$ millisecondes sur un intervalle de temps~$[0,t_{\text{fina}l}]$, et renvoyant la valeur de la moyenne et de l'écart-type des valeurs de ce tableau. On privilégiera l'utilisation de la fonction \texttt{trapezes}.
\end{ques}

\subsection{Validation des mesures}
On sait que chaque moteur de l'imprimante, soumis à un signal d'entrée $e$, doit fournir un signal de sortie $s$ lié au signal d'entrée par l'équation différentielle linéaire du premier ordre :
\begin{equation}
\frac{\mathrm{d} s}{\mathrm{d}t} = -\frac{k}{10}(s-e),
\label{ed}
\end{equation}
où $k$ est un nombre strictement positif. Le bon fonctionnement d'un tel moteur est testé en vérifiant l'adéquation des valeurs mesurées avec des valeurs issues d'une résolution numérique de cette équation.

\begin{ques}
On considère le signal d'entrée $e(t) = \sin(\omega\,t)$ sur l'intervalle de temps $[0,1]$ (en secondes), où~$\omega$ est une constante numérique stockée dans une variable~\texttt{omega}. 

\'Ecrire une fonction \texttt{simulation(...)} renvoyant un tableau de valeurs approchées, calculées à l'aide de la méthode d'Euler, d'une solution de l'équation \eqref{ed} aux instants $[0,0.002,0.004,\dots,1]$ (en secondes). Le candidat sera amené à choisir un ou des arguments pour la fonction programmée, et commentera le choix effectué.

\textbf{On rappelle que tout appel à une fonction de résolution approchée d'équations différentielles  préprogrammée est exclu.}
\end{ques}

On suppose désormais qu'à l'instant $t=0$, la valeur $s(0)$ du signal de sortie doit être nul. Le signal de sortie réel, acquis par le capteur, est stocké dans une variable \texttt{mesures} de type tableau, contenant la valeur du signal aux instants $[0,0.002,0.004,\dots,1]$ (en secondes). L'imprimante doit avoir un comportement conforme à la solution de l'équation \eqref{ed} avec $s(0)=0$. En raison de différences de fabrication, $k$ peut prendre les valeurs $0.5$, $1.1$ ou $2$, et seulement celles-ci. Le comportement du moteur est considéré valide si l'écart maximal entre valeurs observées et valeurs simulées est inférieur à un certain flottant positif \texttt{eps}, pour au moins une valeur de $k$ parmi $0.5$, $1.1$ ou $2$.

\begin{ques}
\'Ecrire une fonction \texttt{validation(mesures,eps)} prenant en arguments un tableau \texttt{mesures} de valeurs observées du signal de sortie et un flottant positif \texttt{eps}, et renvoyant un booléen indiquant si le comportement du moteur peut être validé au seuil \texttt{eps}.
\end{ques}

\section{Bases de données}\label{SQL}
On suppose que les résultats des analyses précédentes ont été stockés dans une base de données, constituée de deux tables, dont on donne une représentation simplifiée dans l'annexe de la page~\pageref{tables}.

Après son assemblage et avant les différents tests de validation, un numéro de série unique est attribué à chaque imprimante. \`A la fin des tests de chaque imprimante, les résultats d'analyse ainsi que le nom du fichier contenant l'ensemble des mesures réalisées sur l'imprimante sont stockés dans la table \texttt{testfin}. Lorsqu'une imprimante satisfait les critères de validation, elle est enregistrée dans la table \texttt{production} avec son numéro de série, la date et l'heure de sortie de production, ainsi que le nom du modèle.

\begin{ques}
\'Ecrire une requête \texttt{SQL} permettant d'obtenir les numéros de série des imprimantes ayant une valeur de \texttt{Imoy} comprise strictement entre $0.4$ et $0.5$.
\end{ques}

\begin{ques}
\'Ecrire une requête \texttt{SQL} permettant d'obtenir les numéros de série, les valeurs de \texttt{Imoy} et \texttt{Iec} ainsi que le modèle des imprimantes ayant été validées et dont la valeur de \texttt{Imoy} est comprise strictement entre 0.4 et 0.5.
\end{ques}

\begin{ques}
\'Ecrire une requête \texttt{SQL} permettant d'obtenir les modèles des imprimantes validées et pour chacun de ces modèles, le nombre d'imprimantes validées.
\end{ques}

\begin{ques}
\'Ecrire une requête \texttt{SQL} renvoyant le numéro de série et le nom du fichier de mesures des imprimantes qui n'ont pas été validées en sortie de production. Commenter très brièvement l'intérêt de cette information pour cette ligne de production.
\end{ques}


\section{Annexe}
\subsection{Base de données}\label{tables}
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline
\multicolumn{5}{|c|}{testfin} \\ \hline
\textbf{nSerie} & dateTest & Imoy & Iec & fichierMes \\ \hline
230-588ZX2547 & 2012-04-22 14-25-45 & 0.45 & 0.11 & mesure31025.csv \\
230-588ZX2548 & 2012-04-22 14-26-57 & 0.43 & 0.12 & mesure41026.csv \\
$\vdots$ & $\vdots$ & $\vdots$ & $\vdots$ & $\vdots$ \\ \hline
\end{tabular}

\vspace*{1ex}
\begin{tabular}{|c|c|c|c|}
\hline
\multicolumn{4}{|c|}{production} \\ \hline
\textbf{Num} & nSerie & dateProd & modele \\ \hline
20 & 230-588ZX2547 & 2012-04-22 15-52-12 & JETDESK-1050 \\
21 & 230-588ZX2549 & 2012-04-22 15-53-24 & JETDESK-3050 \\
$\vdots$ & $\vdots$ & $\vdots$ & $\vdots$\\ \hline
\end{tabular}
\end{center}


  
\end{exercice}





%



\begin{exercice}





\illg{clouseau}{La loi de Newton sur le refroidissement dit que la vitesse de refroidissement d'un objet est
proportionnelle à l'écart de température entre l'objet et le milieu ambiant.\\
L'inspecteur \textsc{Clouseau} arrive sur les lieux d'un meurtre à 9h00. Il commence par prendre la
température de la victime~:~30\degre C. Une heure plus tard, la température du corps est tombée à
29\degre C. Sachant que la température normale du corps d'une personne vivante en bonne santé est de
37\degre C, que la victime était syldave, qu'elle  aimait  les films de gladiateurs et se trouvait
dans une pièce maintenue à 0\degre C, estimez l'heure du décès de la victime et la couleur de ses yeux.}{0.15}{0.8}

On pourra avoir une idée graphique avec nos petites fonctions Python s'inspirant
de la méthode d'Euler.

\end{exercice}



\begin{exercice}




\illt{Dans la forêt syldave, des détritus naturels (~feuilles, branches, animaux morts, cadavres
d'espions, etc.) tombent sur le sol et s'y décomposent. La quantité $Q(t)$ exprimée en $g\cdot m^{-2}$ de
détritus jonchant le sol varie avec le temps $t$. On suppose que de nouveaux détritus tombent au sol
à un taux constant de $200\,g\cdot m^{-2}$ par année et que les détritus accumulés au sol se décomposent
au taux de $50\%$ de la quantité de détritus jonchant le sol.
\begin{enumerate}
\item On note $f(t)=Q(t)-200$. Déterminez une équation différentielle vérifiée par $f$.
\item En déduire l'expression générale de $Q(t)$.
\item Exprimer $Q(t)$ sachant qu'au temps $t=0$, on comptait $50\,g\cdot~m^{-2}$
\end{enumerate}}{paysageSyldave}{0.65}{0.32}


\end{exercice}



%


\begin{exercice}





  \illg{gateway}{
Si  vous allez  vous promener  dans la  ville de  Saint-Louis  dans le
  Missouri aux
  États-Unis, vous  pourrez y admirer  la célèbre \textit{Gateway  Arch to
  the   West}   conçue  en   1947   par   l'architecte  finlandais   Ero
  \textsc{Saarinen}   et   l'ingénieur   Hannskarl  \textsc{Bandel}   et
  dont la construction s'acheva en 1965.  Cette
  arche a la forme d'une \textit{chaînette} pondérée
  d'équation  $$y=212-21\ch(0,033x)$$  où $x$  et  $y$  sont mesurés  en
  mètres  et  rend  hommage  aux  pionniers  partis  à  la  conquête  de
  l'Ouest.\\
Pouvez-vous déterminer  la hauteur de  l'arche et la distance  entre les
deux pieds~?
On rappelle que $\ch x = \frac{\E^x+\E^{-x}}{2}$.}{0.3}{0.65}

\end{exercice}

%

\begin{exercice}




 L'équation   de   la
 hauteur $h$ par rapport au  sol d'un fil électrique suspendu entre deux
 poteaux       s'obtient       en       résolvant       l'\textit{équation
 différentielle}:

    \illt{ \[h''(x)=k \sqrt{1+\left(h'(x)\right)^2}\]
où $k$
 est un paramètre qui dépend de la  densité et de la tension du fil et $x$
est mesuré en  mètres horizontalement à partir d'une  origine située sur
le sol en-dessous du point où la hauteur du fil est la plus faible.

\begin{enumerate}
\item  Vérifiez que  $h~:~x\mapsto  \frac{1}{k}\ch(kx)$ satisfait  cette
      équation différentielle.
\item Quelle  est la hauteur  minimale du fil  si le paramètre  $k$ vaut
      $0,05$~?
\item Quelle  est la  hauteur des poteaux  (de même hauteur)  s'ils sont
      distants de 30~m et que le paramètre $k$ vaut 0,05~?
\end{enumerate}}{cable}{0.65}{0.33}
\end{exercice}




%
\begin{exercice}[DL]
\'{E}crire la formule de Taylor-Mac Laurin et la formule de Taylor-Young 
\`{a} l'ordre $n$  pour les fonctions suivantes et  dans un voisinage de
0:

\begin{multicols}{3}

\begin{enumerate}
\item $f(x)=\text{e}^{x}$

\item $f(x)=\dfrac{1}{1+x}$

\item $f(x)=\dfrac{1}{1-x}$

\item $f(x)=\ln (1+x)$

\item $f(x)=\sin \left( x\right) $

\item $f(x)=\sqrt{1+x},$ $n=4$
\end{enumerate}

\end{multicols}

\end{exercice}


%
\begin{exercice}[DL]
Donner un $DL_{3}(0)$ pour les fonctions suivantes :
\begin{multicols}{4}
\begin{enumerate}
\item $f_{1}(x)=\ln (1+x)$
\item $g(x)=f_{1}(x^{2})$

\item $f_{2}(x)=\text{e}^{x}$

\item $h(x)=f_{2}(-2x)$

\item $f_{3}(x)=f_{1}(x)f_{2}(x)$

\item $f_{4}(x)=f_{1}(x)+f_{2}(x)$
\end{enumerate}
\end{multicols}
\end{exercice}


%

\begin{exercice}[Fonctions homographiques]
  
Une  fonction  homographique est  une  fonction  de  la forme  $x\mapsto
\dfrac{ax+b}{cx+d}$ avec $a$, $b$, $c$ et $d$ des réels.


Notons $g$ la fonction homographique correspondant au cas $a=0$, $b=d=1$
et $c=-1$.

Quel est l'ensemble de définition de $g$~? 

Calculez $g^{(1)}(x)$,  $g^{(2)}(x)$, $g^{(3)}(x$  et plus
généralement $g^{(n)}(x)$.

Écrivez le plus simplement possible $T_{5,0}(g)(x)$.




Sur un  même graphique, donnez  l'allure des courbes  représentatives de
$g$ et de $T_{5,0}(g)$ avec Haskell.



\medskip


Soit $h(x)=\dfrac{5-3x}{1-x}$. Montrez  qu'il existe deux réels $\alpha$
et $\beta$ tels que $h(x)=\alpha+\dfrac{\beta}{1-x}$.

Déduisez-en $T_{5,0}(h)(x)$ pour $x$ appartenant à un bon intervalle.

\end{exercice}











%

\begin{exercice}[Question ouverte...]
  Les polynômes  de \textsc{Newton}  et de  \textsc{Taylor} se  ressemblent non?
  Dans quelle mesure?
\end{exercice}







%1
\begin{exercice}[DS]
  Déterminez  la dérivée  de  la fonction  $x\mapsto \sqrt{\Bigl(\sqrt{\bigl(  \sqrt{(\sqrt
  {x})}\bigr)}\Bigr)}$.  Vous donnerez  l'expression  finale sous  la  forme la  plus
  simple possible.
\end{exercice}


%2
\begin{exercice}[DS]
Déterminez, en détaillant vos calculs, les polynômes de \textsc{Taylor} à l'ordre 4 au voisinage de 0 de:
\begin{multicols}{3}
\begin{enumerate}
\item $f_1(x)=\ln(1-x)$
\item $f_2(x)=\exp(-x)$
\item $f_3(x)=\cos(2x)$
\end{enumerate}
\end{multicols}
\end{exercice}


\begin{exercice}[DS]
 Déterminez  la dérivée  de  la fonction  $t\mapsto \cos{\Bigl(\cos{\bigl(  \cos{(t)}\bigr)}\Bigr)}$.
\end{exercice}


%exo 5

\begin{exercice}[DS]
  Déterminez  le polynôme  de  \textsc{Newton} passant  par  les points  (1,-5),
  (2,-5), (3,1) et (4,19).
\end{exercice}

%exo 6

\begin{exercice}[DS]

On rappelle le résultat suivant:

{\itshape

Soit I un  intervalle, soit $f$ une  fonction de Ιvers $\bbr$  de classe $\CR^n$
avec  $n\in\bbn$   et  soit  $a$  un   élément  de  I.  Il   existe  alors  un
\textit{unique} polynôme  $T_{n,a}$ de degré au  plus $n$ dont les  dérivées jusqu'à
l'ordre $n$ coïncident en $a$ avec celles de $f$.
De plus:
$
T_{n,a}(x) = f(a) + (x-a)f'(a) + \frac{(x-a)^2}{2!}
f''(a)+\frac{(x-a)^3}{3!}f^{(3)}(a)+\cdots +\frac{(x-a)^n}{n!}f^{(n)}(a)
$
  

}

  \begin{enumerate}
  \item Déterminez le  polynôme de \textsc{Taylor} à l'ordre $n$  de la fonction
    ch à partir du théorème rappelé ci-dessus.
  \item Donnez  directement le polynôme de  \textsc{Taylor} à l'ordre $n$  de la
    fonction $t\mapsto \frac{1}{1-t}$.
  \item Déterminez le développement limité à l'ordre 5 au voisinage de 0 de la fonction $t\mapsto
    \frac{\ch(x^2)}{1-3x}$.
  \end{enumerate}
\end{exercice}







%
\begin{exercice}[Courbes de \textsc{Bézier}]

\ill{Dans  les  années 60,  les  ingénieurs  Pierre  \textsc{Bézier} et  Paul
\textsc{De  Casteljau}   travaillant  respectivement  chez   Renault  et
Citroën, réfléchissent  au moyen de  définir de manière la  plus concise
possible la forme d'une carrosserie.

Le  principe  a été  énoncé  par  \textsc{Bézier}  mais l'algorithme  de
construction  par  son  collègue  de  la marque  aux  chevrons  qui  n'a
d'ailleurs été dévoilé  que bien plus tard, la  loi du secret industriel
ayant primé sur le développement scientifique...


Pour la  petite histoire Pierre \textsc{Bézier}  (diplômé de
l'\textsc{Ensam}  et  de \textsc{Supélec}) fut  à  l'origine des  premières
machines à commandes  numériques et de la \textsc{Cao}  ce qui n'empêcha
pas sa  direction de le mettre  à l'écart: il se  consacra alors presque
exclusivement aux
mathématiques  et à  la  modélisation  des surfaces  et  obtint même  un
doctorat en 1977.}{pierre_bezier}

 Paul
\textsc{De Casteljau} était lui un mathématicien d'origine, ancien élève
de la  Rue d'\textsc{Ulm}, qui a  un temps été  employé par l'industrie automobile. 


Aujourd'hui, les courbes de \textsc{Bézier} sont très utilisées en informatique.

Une Courbe de \textsc{Bézier} est une courbe paramétrique aux extrémités imposées
avec des points de contrôle qui définissent les tangentes à cette courbe
à des instants donnés.


\textbf{Algorithme de Casteljau}



Soit $t$ un  paramètre de l'intervalle $[0,1]$ et  $P_1$, $P_2$ et $P_3$
les trois points de contrôle. 

On    construit     le    point    $M_1$     barycentre    du    système
$\left\{(P_1,1-t),(P_2,t)\right\}$    et   $M_2$   celui    du   système
$\left\{(P_2,1-t),(P_3,t)\right\}$. 

On   construit   ensuite   le   point   $M$,   barycentre   du   système
$\left\{(M_1,1-t),(M_2,t)\right\}$.

Exprimez M comme barycentre des trois points $P_1$, $P_2$ et $P_3$.

Faites la construction à la main avec $t=1/3$ par exemple puis utilisez Haskell.



Nous allons assimiler les points $M_1$, $M_2$ et $M$ à des courbes paramétrées. 

Ainsi  $M_1(t)=(1-t)P_1+tP_2$, $M_2(t)=(1-t)P_2+tP_3$ puis

$$M(t)=(1-t)M_1(t)+tM_2(t)=(1-t)^2P_1+2t(1-t)P_2+t^2P_3$$

Vérifiez que $M'(t)=2(M_2(t)-M_1(t))$. Comment l'interpréter~?

\medskip

\textbf{Avec 4 points de contrôle}


Faites une étude similaire (\og à la main\fg{}) avec 4 points de contrôle.

On pourra  utiliser une représentation similaire aux arbres de probabilité.



Créez une  fonction
\verb+bezier3+ qui trace la courbe de \textsc{Bezier} de degré 3 passant par une
liste de 4 points de contrôle.



Avec les points (0,0),(-1,1),(-1,0),(0,-3)



%\begin{figure}
\begin{center}
  \includegraphics{bezier3}
\end{center}
%\end{figure}



Quel est le rôle des différents  points de contrôle?

\medskip


\textbf{Courbe de \textsc{Bézier} du 3\up{e}  degré avec un nombre quelconque de
  points de contrôle}


Il est   pratique  de travailler avec  des polynômes de  degré trois
pour avoir droit à des points d'inflexion.


Augmenter  le  nombre  de  points  de contrôle  implique  a  priori  une
augmentation du degré de la fonction polynomiale.

Pour remédier à  ce problème, on découpe une  liste quelconque en liste
de listes de 4 points. 

Cependant, cela  est insuffisant pour obtenir un  raccordement de classe
$\CR^1$  (pourquoi~? pourquoi est-ce  important d'avoir  un raccordement
de classe $\CR^1$~?)


Pour assurer la  continuité tout court, il faut que  le premier point d'un
paquet soit le dernier du paquet précédent.

Le dernier \og vecteur  vitesse\fg{} de la liste $[P_1,P_2,P_3,P_4]$ est
$\ve{P_3P_4}$. Il  faut donc que ce  soit le premier  vecteur vitesse du
paquet suivant pour assurer la continuité de la dérivée.

Appelons provisoirement le paquet suivant $[P_1',P_2',P_3',P_4']$. On a 
d'une part $P_1'=P_4$
et           d'autre          part          $\ve{P_3P_4}=\ve{P_1'P_2'}$,
i.e. $P_2'=P_4+\ve{P_3P_4}$.

On a donc $P_3'=P_5$ et $P_4'=P_6$.


Connaissant  \verb+bezier3+,  construire  une fonction qui  trace  une
courbe de \textsc{Bézier} cubique avec un nombre quelconque de points de contrôle
(on prendra un nombre pair de points pour se simplifier la vie).

Testez avec:

\begin{pythoncode}
[(1, 1), (2, 3), (5, 5), (6, 2), (7, 7), (10, 5), (10, 2), (8, 0), (4, 0), (5, 1)]
\end{pythoncode}
  


%\begin{figure}
\begin{center}
  \includegraphics{essaiBez}
\end{center}
%\end{figure}


C'était il y a peu  la Saint-Valentin alors
dessinez un c\oe ur....



%\begin{figure}
\begin{center}
  \includegraphics[height=7cm,width=7cm]{coeur}
\end{center}
%\end{figure}



\medskip

\textbf{B-splines uniformes}

Tout ceci est  très beau mais il  y a un hic: en changeant  un point de
contrôle, on modifie grandement la figure.




On considère $m$ n\oe uds $t_0,t_1n,\cdots,t_m$ de l'intervalle $[0,1]$.




Introduisons une nouvelle fonction~:

$$\mathbf{S}(t)= \sum_{i=0}^{m-n-1} \mathbf{P}_{i} b_{i,n}(t) \mbox{ , } t \in [0,1]$$


 les $P_i$ étant les points de contrôle et les fonctions $b_{j,n}$ étant
 définies récursivement par 

$$  b_{j,0}(t)  =  \left\{\begin{matrix}  1  &  \text{si}  \quad  t_j
    \leqslant t < t_{j+1} \\ 0 & \text{sinon} \end{matrix} \right.$$

et pour $n\geqslant  1$

$$    b_{j,n}(t) = \frac{t - t_j}{t_{j+n} - t_j} b_{j,n-1}(t) + \frac{t_{j+n+1} - t}{t_{j+n+1} - t_{j+1}} b_{j+1,n-1}(t). $$

On ne  considérera par  la suite que  des n\oe  uds équidistants~:~ainsi
on aura $t_k=\frac{k}{m}$.

On  parle  de B-splines  uniformes  et  on  peut simplifier  la  formule
précédente en remarquant également des invariances par translation.


À l'aide des formules précédentes, on  peut prouver que dans le cas de 4
points de contrôle on obtient:

$$S(t)=\frac{1}{6}\Big((1-t)^3P_0+(3t^3-6t^2+4)P_1+(-3t^3+3t^2+3t+1)P_2+t^3P_3\Big)$$


Calculez $S(0)$, $S(1)$ puis $S'(0)$ et $S'(1)$:~que peut-on en conclure~?

Reprenez l'étude faite avec les courbes de \textsc{Bézier} 



\begin{pythoncode}
bspline3([(0,0),(-1,1),(-1,0),(0,-3)])
\end{pythoncode}



%\begin{figure}
\begin{center}
  \includegraphics[height=7cm]{bspline3}
\end{center}
%\end{figure}


Obtiendrez-vous un plus joli c\oe ur?


%\begin{figure}
\begin{center}
  \includegraphics[width=8cm,height=8cm]{coeurSpline}
\end{center}
%\end{figure}



\end{exercice}


\modenormal
%%% Local Variables: 
%%% TeX-master: "PolyAnalyse16py"
%%% End: 
