%Last modified: <iut_integral_num.tex modifié par  Guillaume CONNAN le mercredi 16 janvier 2013 à 11h 52min 45s>
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
\chapterimage{ariane}
\chaptertext{À  travers  la  recherche  des  décimales  de  $\pi$,  nous
  plongerons  au  c\oe  ur  du  processeur pour  comprendre  comment  il
  calcule  et   comment  ne  pas  être  surpris   par  quelques  erreurs
  d'arrondis qui pourraient s'avérer catastrophiques...}

\chapter{Intégration numérique et erreurs d'arrondis}


\section{Approximation   de   $\mathbf{\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.



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


\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:


\begin{lstlisting}[caption={}]
def int_rec_droite(f,a,b,N):
    S = 0
    t = a
    dt = (b-a)/N
    while t < b:
        S += f(t)*dt
        t += dt
    return S
\end{lstlisting}

É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{lstlisting}
>>> from math import sqrt,pi
>>> def f(t):
    return sqrt(1-t**2)
>>> for k in range(8): print(pi-4*int_rec_gauche(f,0,1,10**k))
...
3.141592653589793
0.23707432138101003
0.021175621810747725
0.002037186678767622
0.0002011758049373924
2.003710526476965e-05
2.001187481948108e-06
1.9943306694969465e-07
\end{lstlisting}



\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{lstlisting}
>>> for k in range(8):
       print(pi-4*int_trap(f,0,1,10**k))
\end{lstlisting}

Le dernier résultat ne vous paraît-il bizarre?

     
Refaisons alors les calculs mais  avec une boucle  \lstinline+for+ au
lieu  d'une  boucle  \lstinline+while+  et  cette fois  on  obtient  des
résultats plus cohérents. 




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{lstlisting}
>>> for k in range(7): 
        print(pi-12*(int_trap(f,0,1/2,10**k)-sqrt(3)/8))
      \end{lstlisting}
      



Cela semble mieux fonctionner. Bizarre tout de même cette persistance de
\lstinline+1.44+.

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:

\begin{lstlisting}
>>> for k in range(8):
        print(pi-4*int_simps(f,0,1,10**k))
      \end{lstlisting}
      

Pourquoi pas. Et si nous prenions notre petit arc de cercle :

\begin{lstlisting}
>>> for k in range(8): 
        print(pi-12*(int_simps(f,0,1/2,10**k)-sqrt(3)/8))
      \end{lstlisting}
      


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{lstlisting}
>>> from sys import*
>>> float_info.epsilon
2.220446049250313e-16
\end{lstlisting}

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{lstlisting}
>>> from sys import*
>>> e = float_info.epsilon               >>> 0+e/2 == 0
>>> e                                    False
2.220446049250313e-16                    >>> e*1e16
>>> e/2                                  2.220446049250313
1.1102230246251565e-16                   >>> -1-e/2
>>> 1+e/2                                -1.0
1.0                                      >>> 1-e/2
>>> 1+e/2 == 1                           0.9999999999999999
True
\end{lstlisting}

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{lstlisting}
def pseudo_egal_float(a, b):
    return abs(a - b) <= (float_info.epsilon * min(abs(a), abs(b)))
\end{lstlisting}



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

\begin{lstlisting}
>>> pseudo_egal_float(0.1 + 0.1 + 0.1 , 0.3)
True
>>> 0.1 + 0.1 + 0.1 == 0.3                   >>> 1+float_info.epsilon-1
False                                        2.220446049250313e-16
>>> 3 * 0.1 == 0.3                           >>> 10+10*float_info.epsilon-10
False                                        1.7763568394002505e-15
>>> 4 * 0.1 == 0.4                           >>> x = 0.1
True                                         >>> 3*x-0.3
>>> 10+float_info.epsilon-10                 5.551115123125783e-17
0.0                                          >>> 4*x-0.4
                                             0.0
\end{lstlisting}


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; posons 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 alors  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{lstlisting}
>>> for k in range(13): 
        print(pi-12*(int_simps(f,0,1/2,2**k)-sqrt(3)/8))
      \end{lstlisting}
      


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{lstlisting}
>>> f(1e-8)
1.0
\end{lstlisting}



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{lstlisting}
>>> 1+0.1+0.1+0.1-0.3<=1
False
\end{lstlisting}


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{lstlisting}
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)))
\end{lstlisting}


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}.


\subsection{Quelques mots sur la manipulations des flottants}

En 1985, une  norme de représentations des nombres  a été proposée afin,
entre autre, de  permettre de faire des programmes  portables: il s'agit
de la norme IEEE-754 (Standard for Binary Floating-Point Arithmetic).

En simple précision, un nombre est représenté sur 32 bits (en fait 33...):
\begin{itemize}
\item le premier donne le signe;
\item un 1 implicite puis 23 bits pour la partie fractionnaire;
\item 8 bits pour l'exposant.
\end{itemize}


Les  24 bits  qui ne  sont ni  le signe,  ni l'exposant  représentent la
mantisse qui appartient à $[1,2[$.

Notons $s_x$ le bit de  signe, $e_x$ les bits d'exposant, $m_x=1+f_x$ la
mantisse avec $f_x$ la partie fractionnaire.

Il y a deux zéros, $-0$ et $+0$ car tout nombre est signé. Par exemple:

\begin{itemize}
\item $s_{-0}=1$;
\item $ b_{-0}=00000000$
\item $f_{-0}=00000000000000000000000$
\end{itemize}

Il   y   a   aussi    deux   infinis.   Par   exemple   $s_{-\infty}=1$,
$e_{-\infty}=11111111$ et $f_{-\infty}=0...$.

Il y  a enfin  le \textit{Not a  Number}, noté  NaN, qui est  codé comme
$+\infty$ mais avec une partie fractionnaire non nulle.


\medskip


La  majorité  des  résultats   troublants  observés  vient  des  erreurs
d'arrondis,  les  deux  principales  étant l'\textit{élimination}  et
l'\textit{absorption}.

La  première intervient  lors de  la soustraction  de deux  nombres très
proches. Par exemple:

\begin{verbatim}
 1.10010010000111111011011
-1.10010010000110000000000
--------------------------
=0.00000000000111111011011
\end{verbatim}

La mantisse est ensuite renormalisée pour devenir:

\begin{verbatim}
 1.11111011011000000000000
\end{verbatim}

Les zéros ajoutés à droite sont faux.

\smallskip

L'absorption  intervient lorsqu'on  additionne deux  nombres  d'ordre de
grandeur très différents: les informations concernant le plus petit sont
perdues.


\begin{verbatim}
 1.100100100001111100 01011
+                   1.00000001111
---------------------------------
=1.10010010000111110100101101111
\end{verbatim}

Mais après normalisation, on a:

\begin{verbatim}
1.10010010000111110100101
\end{verbatim}

Voici     un      exemple     classique:     le      calcul     de     $
\displaystyle\sum_{k=1}^{2^{15}}\frac{1}{k^2}$ dans un sens puis dans l'autre.


\medskip

\begin{multicols}{2}
\begin{lstlisting}
def som_croissante(N):
    S=0
    for k in range(1,N+1):
        S+=1/k**2
    return S
  \end{lstlisting}
  
\begin{lstlisting}
def som_decroissante(N):
    S=0
    for k in range(N,0,-1):
        S+=1/k**2
    return S
\end{lstlisting}
\end{multicols}

Alors:

\begin{lstlisting}
>>> som_croissante(2**15)
1.6449035497357558
>>> som_decroissante(2**15)
1.6449035497357580
\end{lstlisting}

En fait, le résultat correct est:

\begin{lstlisting}
sage: k=var('k')
sage: s=sum((1/k**2),k,1,2**15)
sage: s.n(digits=20)
1.6449035497357579868
\end{lstlisting}





\medskip

Les conséquences de  telles erreurs non prises en  compte peuvent être plus
graves qu'une séance de TP ratée.

Le 4 juin 1996  par exemple, la fusée Ariane 5 a  explosé en vol, trente
secondes  après son  décollage.  Après  enquête, il  s'est  avéré que  la
vitesse horizontale  de la fusée par  rapport au sol  était calculée sur
des  flottants 64 bits  puis convertie  en entier  signé 16  bits. Cette
méthode  avait été appliquée  sur Ariane  4 avec  succès car  sa vitesse
était plus  faible donc tenait  sur un entier  16 bits, mais  ce n'était
plus le cas  pour Ariane 5...\footnote{Voir le rapport  de la commission
  présidée  par  J.L.   \textsc{Lions}  disponible  à  cette  adresse:\\
  \href{http://www.ima.umn.edu/~arnold/disasters/ariane5rep.html}{http://www.ima.umn.edu/~arnold/disasters/ariane5rep.html}
  ainsi que la page 22 d'un article de William \textsc{Kahane} : \href{http://www.eecs.berkeley.edu/~wkahan/JAVAhurt.pdf}{http://www.eecs.berkeley.edu/~wkahan/JAVAhurt.pdf}}

Un problème d'élimination dans l'horloge des missiles \textit{Patriot} a
causé également  la mort  de dizaines de  personnes pendant  la première
guerre du Golfe\footnote{Voir \href{http://ta.twi.tudelft.nl/nw/users/vuik/wi211/disasters.html}{http://ta.twi.tudelft.nl/nw/users/vuik/wi211/disasters.html}}.


%\subsection{Utilisation des ressources \textsf{Python}{} }



\begin{remarque}

Les erreurs présentées ici ne sont pas dues à \textsf{Python}{}et on les retrouve sur
d'autres langages:

\begin{itemize}

\item avec CAML:

\begin{Verbatim}
# 0.3-.3.*.0.1;;
- : float = -5.5511151231257827e-17
\end{Verbatim}

\item avec SCILAB:

  \begin{Verbatim}
-->0.3-3*0.1
 ans = - 5.551D-17
  \end{Verbatim}

\item avec giac/XCAS:

  \begin{Verbatim}
1>> 0.3-3*0.1
1.42108547152e-14
// Time 0
  \end{Verbatim}

\item avec MAXIMA:

  \begin{Verbatim}
(%i1) 0.3-3*0.1;
(%o1) - 5.5511151231257827E-17
  \end{Verbatim}


\end{itemize}

\end{remarque}


\begin{recherche}
On   considère   la  suite   $u$   définie   par  $u_{n+1}=4u_n-1,\   \
u_0=\frac{1}{3}$.

Calculez à la main les premiers termes de la suite.

Utilisez ensuite Python pour obtenir une approximation des 100 premiers termes.

Faites de même avec la suite $v$ définie par $v_{n+1}=3v_n-1,\ \ v_0=0,5$.
\end{recherche}

\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.

Programmons les trois méthodes sus-mentionnées :

\begin{lstlisting}
class Rectangles(Integration):
    def quadrature(self):
        a, b, n = self.a, self.b, self.n
        h = (b-a) / n
        return [(h, a + i*h) for i in range(n)]

class Trapezes(Integration):
    def quadrature(self):
        a, b, n = self.a, self.b, self.n
        h = (b-a) / n
        return [(h, a + (i+0.5)*h) for i in range(n)]

class Simpson(Integration):
    def coeff(self, i, N):
        if i == 0 or i == N:
            return 1
        elif i % 2 == 1:
            return 4
        else:
            return 2

    def quadrature(self):
        a, b, n = self.a, self.b, self.n
        if n % 2 != 1:
            n += 1
        h = (b-a) * 0.5 / n
        return [(h/3*self.coeff(i, 2*n), a + i*h) for i in range(2*n+1)]
\end{lstlisting}


Pour diminuer le nombre de points à tracer, et donc le temps de calcul, on
pourrait utiliser des courbes de \textsc{Bézier} avec comme points de contrôle les
vecteurs vitesse.

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 $\mathbf{\pi}$ et les arctangentes}

\subsection{Le nombre $\mathbf{\pi}$ et \textsc{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{lstlisting}
>>> from decimal import*
# on regle la precision a 30 chiffres
>>> getcontext().prec = 30
# on rentre les nombres entre apostrophes, comme des chaines
>>> deux=Decimal('2')
>>> deux.sqrt()
Decimal('1.41421356237309504880168872421')
>>> getcontext().prec = 50
>>> deux.sqrt()
Decimal('1.4142135623730950488016887242096980785696718753769')
\end{lstlisting}




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.



\begin{lstlisting}
def greg(a,N):
    ad = Decimal('1')/Decimal(str(a))
    nd = ad
    dd = Decimal('1')
    s = nd/dd
    for k in range(1,N):
        nd *= Decimal('-1')*ad*ad
        dd += Decimal('2')
        s += nd/dd
    return s
\end{lstlisting}

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


\begin{lstlisting}[caption={}]
sage: pi.n(digits=1000)
\end{lstlisting}

On la copie dans \textsf{Python}{} et on compare:

\begin{lstlisting}
>>> PI-pi_machin(1000)
\end{lstlisting}

En une seconde et demie, on a mille bonnes décimales de $\pi$...
% 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!






%%% Local Variables: 
%%% TeX-master: "poly_info2_asymptotic"
%%% End: 
