TP IPT 5 : problèmes de floataison

Kerviel ? Un amateur !

M. X a récemment été à sa banque (Chaotic Bank Society), pour connaître les nouvelles offres proposées aux
meilleurs clients. Son banquier lui propose l'offre suivante: « vous déposez tout d'abord e - 1 euros sur votre compte
(où e = 2.7182818... est la base du logarithme népérien). La première année, nous prenons 1 euro sur votre compte de
frais de gestion. Par contre, la deuxième année est plus intéressante pour vous, car nous multiplions votre capital restant
par 2 et prenons 1 euro de frais de gestion. La troisième année est encore plus intéressante, car nous multiplions votre
capital par 3 et prenons 1 euro de frais de gestion. Et ainsi de suite : la n-ième année, nous multiplions votre capital par n
et prenons 1 euro de frais de gestion. Intéressant, non? » Pour prendre sa décision, M. X décide de demander l'aide
d'un informaticien et se retourne vers vous. (Exemple dû à JM Muller)

Vous ferez une belle simulation numérique du problème agrémentée de beaux graphiques :

Python
from math import exp
import matplotlib.pyplot as plt
 
plt.plot(les abscisses, les ordonnées)
plt.show() # je regarde
plt.clf() # j'efface

Une petite étude mathématique d'une certaine suite s'avère peut-être nécessaire pour éviter d'être pris pour un gogo.

Les nombres plus petits que 1 en base 2 ou 16

Déterminez deux fonctions qui renvoient l'écriture en base 2 puis en base 16 d'un nombre positif plus petit que 1.

Python
In [1]: decToBin(0.375,10)
Out[1]: '0b0.0110000000'
In [2]: decToHex(0.90625,3)
Out[2]: '0x0.e80'

On rappelle que la fonction chr renvoie le caractère correspondant au code ASCII donné en argument:

Python
In [3]: chr(97)
Out[3]: 'a'

On pourra commencer par donner une fonction qui renvoie une chaîne de caractère correspondant à la représentation en base 16 d'un entier:

Python
In [4]: hstr(1000)
Out[4]: '3e8'

Quand on ne sait pas compter...

Le 25 février 1991, un missile Patriot a raté son Scud qui a tué 28 soldats américains et blessé cent autres parce que son horloge comptait le temps en 1/10e de seconde.
Représentez 1/10 avec une mantisse de 24 bits et tronquez le résultat.

Calculez l'erreur en seconde puis après 100 heures d'utilisation. Sachant qu'un
Scud vole à ${1676}{m.s^{-1}}$, quelle est environ l'erreur commise en mètres?

Et si on avait arrondi au plus proche?

Calcul de $\sqrt{2}$

Il y a presque 4000 ans, les petits Babyloniens calculaient la racine carrée de
2 avec 6 bonnes décimales en utilisant un algorithme bien connu.

$$r_{n+1}= f(r_n)=\frac{r_n + \frac{2}{r_n}}{2},\qquad r_0=1$$

Est-ce que la suite des valeurs calculées par la boucle
converge vers $\sqrt{2}$?

Démontrez que pour tout entier naturel non nul $n$, on a $r_n \geqslant \sqrt{2}
$ puis que la suite est décroissante.

quelle est la vitesse de convergence? Ici, cela se
traduit par « combien de décimales obtient-t-on en plus à chaque
itération?».

Calculez $(r_{n+1} - \sqrt{2} )$ en fonction de $(r_n- \sqrt{2} )$ et essayez
de répondre à la question précédente à l'aide de ce résultat.

Soit $(r_n)$ une suite convergeant vers $\ell$. S'il existe un entier $k>0$ tel que:

$$
\displaystyle\lim_{n\to +\infty}\frac{|r_{n+1}-\ell|}{|r_n-\ell|^k}=C
$$

avec $C\neq 0$ et $C\neq +\infty$ alors on dit que $(r_n)$ est d'ordre $k$ et
que C est la constante asymptotique d'erreur.

Déterminez par exemple l'ordre et la constante asymptotique d'erreur de la suite
de Babylone donnant une approximation de $\sqrt{2}$.

Si l'on connaît l'ordre d'une suite, on va pouvoir prévoir l'évolution de la
précision atteinte au fur et à mesure de l'avance de la boucle.

Ainsi, pour une suite d'ordre $k$ avec comme constante asymptotique d'erreur C,
on a pour $n$ suffisamment grand:

$$|r_{n+1}-\ell|\approx C|r_n-\ell|^k$$

En posant $d_n=-\log_{10}|r_n-\ell|$ et $D=-\log_{10}C$ on obtient:

$$d_{n+1}=...$$

Que peut-on en conclure?

On pourra étudier également un exercice du poly de M. Sauvageot.

Approximation de $\pi$

Donnez une approximation de $\pi$ à l'aide d'une approximation de calcul d'intégrale.

Python
In [30]: cintegre(0,1,100000000) - pi
Out[30]: 1.0476064460362977e-12
 
In [31]: cintegre(0,1,1000000000) - pi
Out[31]: -1.0480505352461478e-13

Équation du second degré

Vous savez résoudre une équation du type $ax^2+bx+c=0$ avec
$a\neq 0$...

Les racines, si elles existent, sont données par une formule bien connue
dépendant de $a$, $b$ et $c$:

$$r=\frac{-b\pm \sqrt{b^2-4ac}}{2a}=\frac{-b\pm \sqrt{Δ}}{2a} \quad
\text{avec}\quad Δ = b^2-4ac$$

Où peuvent se cacher d'éventuelles annulations catastrophiques? Étudiez ces cas
avec attention, voyez si vous pouvez éviter les éliminations catastrophiques en
réécrivant les formules un peu dans l'esprit de la « levée
d'indétermination ».

  1. Que se passe-t-il lorsque $b^2 \gg |4ac|$? En quoi la formule $\frac{-
    \left(b+\operatorname{signe}(b)\sqrt{Δ}\right)}{2a} $ peut aider?
  2. Que se passe-t-il lorsque $b^2\approx 4ac$? Peut-on y remédier? Que
    peut-on dire de $Δ$ par rapport à $b^2$? Y a-t-il élimination catastrophique?
  3. Que se passe-t-il dans le cas de l'équation $10^{200}x^2-3\times
    10^{200}x+2\times 10^{200}=0$?
  4. Et dans le cas de $10^{-200}x^2-3x+2\times 10^{200}=0$?
  5. Moralité?

Une proposition de correction

Python
from math import exp, sqrt, pi
import matplotlib.pyplot as plt
import cython
 
def decToBin(x,n):
    y = x
    s = "0b0."
    for k in range(n):
        y *= 2
        b = int(y)
        s += str(b)
        y -= b
    return s
         
 
def hstr(d):
    if d < 10:
        return str(d)
    elif d < 16:
        return chr(97 + d - 10)
    else:
        return (hstr(d // 16) + hstr(d % 16))
 
def decToHex(x,n):
    y = x
    s = "0x0."
    for k in range(n):
        y *= 16
        b = int(y)
        s += hstr(b)
        y -= b
    return s
    
def chaotic(n,c0):
    cpt = c0
    for i in range(1,n + 1):
        print(cpt)
        cpt = i*cpt - 1
    return cpt
 
# Avec numpy :
# import numpy as np
# listeC0 = np.linspace(1.66,1.76,1000)
# p = plt.plot(listeC0,chaotic(listeC0,50))
# plt.show()
# ou, sans numpy :
# listeC0 = [1.66 + k*0.001 for k in range(101) ]
# listeC50 = list(map(lambda c: chaotic(50,c), listeC0))
# p = plt.plot([eC0,listeC50])
    
def heron(x,n):
    if n == 0:
        return x
    return heron((x + 2/x)*0.5, n - 1)
 
 
def c(x):
    return 4*sqrt(1 - x**2)
 
def rect(f,a,b,n):
    s = 0
    h = (b - a) / n
    for k in range(n):
        s += h * f(a + (k + 0.5)*h)
    return s
 
# %%cython
# cdef extern from "math.h":
#     double c_sqrt "sqrt"(double)
 
# cdef double cf(double x):    
#    return 4*c_sqrt(1 - x**2)
 
 
# def cintegre(double a, double b, long N):
#     cdef int i
#     cdef double s, dx
#     s  = 0
#     dx = (b - a) / N
#     for i from 0 <= i < N:
#         s += cf(a + (i + 0.5)*dx)
#     return s*dx
 
 
# #In [31]: cintegre(0,1,10**9) - pi
# #Out[31]: -1.0480505352461478e-13

[/]