Chapitre 34 — calculer les termes d'une suite, programmer un algorithme de seuil, résoudre $f(x)=0$ par dichotomie, sécante et Newton, et comprendre pourquoi $0{,}1+0{,}2\neq 0{,}3$ en machine.
Naviguez avec les flèches ← → du clavier, ou via le Sommaire.
for, liste des termesnumpyPython conjecture et guide la démonstration ; seule la démonstration établit le résultat.
Si $u_n=f(n)$, la fonction u(n) renvoie directement $f(n)$. Si $u_{n+1}=g(u_n)$, elle part de $u_0$ et applique $g$ à $n$ reprises, avec une boucle for de $n$ tours.
def u(n) : x = 1 # u_0 for i in range(n) : x = x/2 + 3 # de u_i à u_{i+1} return x >>> u(10) 5.9951171875 >>> u(60) 6.0
On démontre que $u_n=6-\dfrac{5}{2^n}\neq 6$, et pourtant Python affiche 6.0 pour $u_{60}$ : c'est un effet des nombres flottants (voir plus loin).
def fibo(n) : a, b = 0, 1 # F_0, F_1 for i in range(n) : a, b = b, a+b # F_{i+1}, F_{i+2} return a def termes(n) : liste = [1] # u_0 for i in range(n) : liste.append(liste[i]/2 + 3) return liste >>> fibo(100) 354224848179261915075 >>> termes(5) [1, 3.5, 4.75, 5.375, 5.6875, 5.84375]
a, b = b, a+b calcule d'abord les deux valeurs de droite : pas de variable temporaire.termes(n)[k] est $u_k$ ; pour représenter la suite, un nuage de points : plt.plot(x, y, 'o').while.n : u contient toujours $u_n$.u = 1 n = 0 while u <= 2 : # capital placé à 5 % u = 1.05*u n = n+1 print(n, u) 15 2.0789281794113688
Le capital a doublé au bout de $15$ ans : $u_{14}\leqslant 2<u_{15}$.
def seuil(eps) : # u_0 = 1, u_{n+1} = u_n/2 + 3, limite 6 u = 1 n = 0 while abs(u-6) >= eps : u = u/2 + 3 n = n+1 return n >>> seuil(0.01) 9 >>> seuil(1e-6) 23
Si $(u_n)$ est monotone et converge vers $\ell$, la boucle se termine et le rang renvoyé est exactement le $n_0$ de la définition de la limite : $\lvert u_n-\ell\rvert<\varepsilon$ pour tout $n\geqslant n_0$.
Vérification : $\dfrac{5}{2^n}<\varepsilon\iff 2^n>\dfrac{5}{\varepsilon}$ ; pour $\varepsilon=10^{-2}$, $2^8<500<2^9$, d'où $n=9$.
Si la suite ne tend pas vers $\ell$, la boucle peut ne jamais s'arrêter ; seuil(0) tourne indéfiniment. Sans monotonie, un terme ultérieur peut ressortir de la bande.
from math import sqrt u = 0 # u_{n+1} = sqrt(2 + u_n) for n in range(1, 9) : u = sqrt(2+u) print(n, u) 1 1.4142135623730951 2 1.8477590650225735 3 1.9615705608064609 4 1.9903694533443939 5 1.9975909124103448 6 1.9993976373924085 7 1.999849403678289 8 1.9999623505652022
La suite semble croissante, majorée par $2$, de limite $2$. Par récurrence, $0\leqslant u_n\leqslant u_{n+1}\leqslant 2$ : elle converge vers $\ell\in[0,2]$, et la continuité de $x\longmapsto\sqrt{2+x}$ donne $\ell=\sqrt{2+\ell}$, soit $\ell^2-\ell-2=0$, donc $\ell=2$.
h = 0 # H_n = 1 + 1/2 + ... + 1/n for k in range(1, 10**6+1) : h = h + 1/k print(h) 14.392726722864989
Un million de termes donnent $H_n\simeq 14{,}39$ : la suite semble s'essouffler. Pourtant $H_{2n}-H_n\geqslant\dfrac12$, donc $H_n\to+\infty$. Un calcul suggère, il ne prouve pas.
Si $(u_n)$ croît, $(v_n)$ décroît et qu'elles sont adjacentes, $u_n\leqslant\ell\leqslant v_n$ : encadrement d'amplitude $v_n-u_n\to 0$. Ici :
def approx_e(p) : n, fact, u = 1, 1, 2 # u_1 = 2 v = u + 1 # v_1 = 3 while v-u > 10**(-p) : n = n+1 fact = fact*n # n! mis à jour, pas recalculé u = u + 1/fact v = u + 1/(n*fact) return n, u, v >>> approx_e(10) (13, 2.7182818284467594, 2.7182818284591126)
u = 2 # u_0 for n in range(6) : v = 2/u # v_n print(n, v, u) u = (u + 2/u)/2 # u_{n+1} 0 1.0 2 1 1.3333333333333333 1.5 2 1.411764705882353 1.4166666666666665 3 1.41421143847487 1.4142156862745097 4 1.4142135623715002 1.4142135623746899 5 1.4142135623730951 1.414213562373095
$v_n\leqslant\sqrt2\leqslant u_n$ : l'amplitude passe de $1$ à $3\cdot10^{-12}$ en quatre étapes, le nombre de décimales exactes double environ à chaque fois.
L'affichage donne $v_5>u_5$, ce qui est mathématiquement impossible : on a atteint la limite de précision de la machine.
Un int est exact et illimité. Un float est stocké sur $64$ bits sous la forme $\pm m\times 2^{e}$, avec une mantisse $m$ de $53$ chiffres binaires et $-1022\leqslant e\leqslant 1023$ : environ $16$ chiffres significatifs, valeur absolue au plus $1{,}8\cdot10^{308}$.
>>> 0.1 + 0.2 0.30000000000000004 >>> 0.1 + 0.2 == 0.3 False >>> 0.5 + 0.25 == 0.75 True >>> (0.1).as_integer_ratio() (3602879701896397, 36028797018963968)
En base $2$, $\dfrac1{10}=0{,}0001100110011\dots$ a une écriture infinie : Python stocke le flottant le plus proche (une fraction de dénominateur $2^{55}$). $0{,}5$ et $0{,}25$, eux, sont exacts.
>>> 2**100 1267650600228229401496703205376 >>> 2.0**100 1.2676506002282294e+30 >>> 2.0**1100 OverflowError: (34, 'Result too large') >>> 10**16 + 1 - 10**16 1 >>> 1e16 + 1 - 1e16 0.0
int, + - * ** // % sont exacts./ renvoie toujours un flottant ; un seul opérande flottant suffit pour un résultat flottant, arrondi à $53$ bits.Près de $10^{16}$, deux flottants consécutifs sont distants de $2$ : $10^{16}+1$ est arrondi à $10^{16}$. De même $6-\dfrac{5}{2^{60}}$ est arrondi à $6$, d'où u(60) qui affiche 6.0.
s = 0 for i in range(10) : s = s + 0.1 print(s, s == 1) 0.9999999999999999 False >>> abs(0.1 + 0.2 - 0.3) < 1e-9 True
== entre flottants calculés : abs(a-b) < eps ou isclose(a, b) du module math.while x < 1, jamais while x != 1 (boucle infinie en ajoutant $0{,}1$).f(a)*f(m) <= 0, pas f(m) == 0.>>> round(2.675, 2) 2.67 >>> round(2.5), round(3.5) (2, 4) >>> x = 123456.789 >>> f"{x:.3e}", f"{x:.1f}" ('1.235e+05', '123456.8')
round(x, k) arrondit à $k$ décimales, int(x) tronque. Le flottant stocké pour $2{,}675$ vaut $2{,}67499999\dots$, d'où $2{,}67$.
round arrondit à l'entier pair le plus proche (norme IEEE 754). Au-delà de la seizième décimale, les chiffres affichés n'ont aucun sens.
# f(x) = x**3 + x - 1, unique solution c dans [0, 1] def dicho2(f, a, b, p) : n = 0 # nombre d'itérations fa = f(a) while b-a > 10**(-p) : m = (a+b)/2 fm = f(m) if fa*fm <= 0 : b = m else : a, fa = m, fm n = n+1 return (a+b)/2, n >>> dicho2(f, 0, 1, 10) (0.6823278038355056, 34)
Après $n$ itérations, l'intervalle mesure $\dfrac{b-a}{2^n}$ : la boucle s'arrête pour
Environ $3{,}3$ itérations par décimale. Le balayage de pas $h$ coûte, lui, $\dfrac{b-a}{h}$ évaluations.
La droite passant par $\big(x_0,f(x_0)\big)$ et $\big(x_1,f(x_1)\big)$ coupe l'axe des abscisses en
On recommence avec $x_1$ et $x_2$, et ainsi de suite.
def secante(f, x0, x1, eps) : n = 0 while abs(x1-x0) > eps : x2 = x1 - f(x1)*(x1-x0)/(f(x1)-f(x0)) x0, x1 = x1, x2 n = n+1 return x1, n >>> secante(f, 0, 1, 1e-10) (0.6823278038280193, 8)
La tangente en $x_n$ coupe l'axe des abscisses en (si $f'(x_n)\neq 0$)
Si $x_0$ est assez proche de $c$, convergence quadratique : les décimales exactes doublent à chaque itération.
x = 1 # df(x) = 3*x**2 + 1 for n in range(1, 7) : x = x - f(x)/df(x) print(n, x) 1 0.75 2 0.686046511627907 3 0.6823395825973142 4 0.6823278039465127 5 0.6823278038280194 6 0.6823278038280193
>>> import numpy as np >>> t = np.array([1, 2, 3, 4]) >>> t**2 array([ 1, 4, 9, 16]) >>> np.linspace(0, 1, 5) array([0. , 0.25, 0.5 , 0.75, 1. ]) >>> n = np.arange(0, 5) >>> n**2 / 2**n array([0. , 0.5 , 1. , 1.125, 1. ])
Opérations et fonctions (np.sqrt, np.exp…) agissent terme à terme, sans boucle. Pour une courbe : x = np.linspace(0, 1, 200) puis plt.plot(x, f(x)).
Les entiers numpy tiennent sur $64$ bits : np.array([2])**63 renvoie $-9223372036854775808$ sans erreur. Pour de grands entiers, listes et entiers Python.
for de $n$ tourswhile + négation de la conditionabs(a-b) < eps, jamais ==while.n avant de calculer le terme : u et n se décalent d'un rang.0.1 + 0.2 == 0.3 : c'est False.dicho(f,0,1,17) ne se termine jamais.while : le nombre de tours n'est pas connu à l'avance.0.1 + 0.2 == 0.3 ?False : $0{,}1$ n'a pas d'écriture binaire finie.