Chapitre 35 — approcher une intégrale, suivre une solution d'équation différentielle point par point, et simuler le hasard pour confronter fréquences et probabilités.
Naviguez avec les flèches ← → du clavier, ou via le Sommaire.
random, Bernoulli, loi binomialeQuand la formule manque, on calcule : une somme approche une intégrale, une ligne brisée approche une courbe, une fréquence approche une probabilité.
$h=\dfrac{b-a}{n}$, $x_k=a+kh$. Rectangles à gauche et à droite :
from math import * def rect(f,a,b,n) : h=(b-a)/n u,v=0,0 x=a for i in range(n) : u=u+h*f(x) x=x+h v=v+h*f(x) return (u,v) def f(x) : return log(1+x) >>> rect(f,0,1,10) (0.3512205777177568, 0.4205352957737513) >>> rect(f,0,1,100) (0.38282445857472913, 0.3897559303803286) >>> rect(f,0,1,1000) (0.3859477458629464, 0.38664089304350635) >>> 2*log(2)-1 0.3862943611198906
Si $f$ est croissante sur $[a,b]$, $S_n\leqslant\displaystyle\int_a^b f(x)\,dx\leqslant S_n'$ (renversé si $f$ décroît), et l'amplitude vaut
Pour $\displaystyle\int_0^1\ln(1+x)\,dx=2\ln2-1$ à $10^{-3}$ près : $n\geqslant\ln2\times10^3\simeq693{,}1$, soit $n=694$. Chaque multiplication de $n$ par $10$ n'apporte qu'une décimale : la méthode est lente.
>>> rect(f,0,1,694) (0.3857948890324332, 0.3867936601859778)
Hauteur prise au milieu de chaque intervalle, ou trapèze dont deux sommets sont sur la courbe : $T_n=\dfrac{S_n+S_n'}{2}$.
def milieu(f,a,b,n) : h=(b-a)/n s=0 for i in range(n) : s=s+f(a+(i+0.5)*h) return h*s def trapezes(f,a,b,n) : h=(b-a)/n s=(f(a)+f(b))/2 for i in range(1,n) : s=s+f(a+i*h) return h*s >>> milieu(f,0,1,10) 0.38650248251865776 >>> trapezes(f,0,1,10) 0.38587793674575405
I=2*log(2)-1 def err(x) : return f"{abs(x-I):.2e}" for n in [10,20,40,80,160] : u,v=rect(f,0,1,n) m=milieu(f,0,1,n) t=trapezes(f,0,1,n) print(n,err(u),err(m),err(t)) # Sortie : n, rectangles, milieu, trapèzes 10 3.51e-02 2.08e-04 4.16e-04 20 1.74e-02 5.21e-05 1.04e-04 40 8.69e-03 1.30e-05 2.60e-05 80 4.34e-03 3.26e-06 6.51e-06 160 2.17e-03 8.14e-07 1.63e-06
Quand $n$ double, l'erreur des rectangles est divisée par $2$ (ordre $\frac1n$), celle du point milieu et des trapèzes par $4$ (ordre $\frac1{n^2}$).
Point milieu par excès, trapèzes par défaut ici ($\ln$ est concave). Aucune des deux ne donne d'encadrement garanti.
Pour $y'=f(x,y)$, $y(x_0)=y_0$ et $h$ petit, la courbe est proche de sa tangente sur $[x,x+h]$ :
Le schéma d'Euler de pas $h$ est la suite de points $(x_k,y_k)$ :
$y_k$ approche $y(x_k)$ ; la ligne brisée approche la courbe de la solution.
def euler_affine(a,b,x0,y0,h,n) : X=[x0] Y=[y0] x,y=x0,y0 for k in range(n) : y=y+h*(a*y+b) x=x+h X.append(x) Y.append(y) return X,Y X,Y=euler_affine(1,0,0,1,0.1,10) # y'=y, y(0)=1 for k in range(0,11,2) : print(round(X[k],1), round(Y[k],4), round(exp(X[k]),4), round(abs(Y[k]-exp(X[k])),4)) # Sortie : x, y_k, exp(x), écart 0 1 1.0 0.0 0.2 1.21 1.2214 0.0114 0.4 1.4641 1.4918 0.0277 0.6 1.7716 1.8221 0.0506 0.8 2.1436 2.2255 0.082 1.0 2.5937 2.7183 0.1245
L'erreur croît le long de l'intervalle : chaque pas part d'un point déjà faux, et la tangente, en-dessous d'une courbe convexe, sous-estime la solution.
for n in [10,100,1000,10000] : X,Y=euler_affine(1,0,0,1,1/n,n) print(n,1/n,Y[-1],abs(Y[-1]-e)) # Sortie : 10 0.1 2.5937424601 0.124539368359045 100 0.01 2.704813829421526 0.013467999037519274 1000 0.001 2.716923932235896 0.0013578962231490799 10000 0.0001 2.7181459268252266 0.00013590163381849152
À abscisse fixée, l'erreur est proportionnelle au pas $h$ : diviser $h$ par $10$ divise l'erreur par $10$, pour dix fois plus de calculs.
Pour $y'=y$ : $y_{k+1}=(1+h)y_k$, donc avec $h=\frac1n$, $y_n=\left(1+\frac1n\right)^n\longrightarrow e$.
def euler(f,x0,y0,h,n) : X=[x0] Y=[y0] x,y=x0,y0 for k in range(n) : y=y+h*f(x,y) x=x+h X.append(x) Y.append(y) return X,Y def f(x,y) : return -2*x*y X,Y=euler(f,0,1,0.1,20) for k in range(0,21,4) : print(round(X[k],1), round(Y[k],4), round(exp(-X[k]**2),4)) # Sortie : x, y_k, exp(-x²) 0 1 1.0 0.4 0.8844 0.8521 0.8 0.5542 0.5273 1.2 0.2382 0.2369 1.6 0.0675 0.0773 2.0 0.012 0.0183
Solution exacte $x\mapsto e^{-x^2}$. En $x=1$, l'erreur passe de $1{,}4\times10^{-2}$ à $1{,}2\times10^{-3}$ puis $1{,}2\times10^{-4}$ pour $h=0{,}1$, $0{,}01$, $0{,}001$.
randomrandom() : flottant de $[0,1[$, loi uniforme ; uniform(a,b) : flottant de $[a,b]$randint(a,b) : entier de $\llbracket a,b\rrbracket$, bornes incluseschoice(L) : un élément de L ; shuffle(L) : mélange sur place, ne renvoie rienseed(g) : fixe la graine, donc toute la suite de tirages>>> import random >>> random.seed(2026) >>> random.random() 0.11911988496396309 >>> random.randint(1,6) 5 >>> random.choice(["pile","face"]) 'pile' >>> L=[1,2,3,4,5] >>> random.shuffle(L) >>> L [1, 5, 3, 4, 2] >>> random.uniform(2,5) 4.253077874475854
random() suit la loi uniforme sur $[0,1[$, donc l'événement random()<p a pour probabilité $p$ : c'est le succès. Pièce truquée à $0{,}7$ : bernoulli(0.7).
def bernoulli(p) : if random.random()<p : return 1 else : return 0 def de() : return random.randint(1,6) >>> random.seed(1) >>> [bernoulli(0.7) for i in range(10)] [1, 0, 0, 1, 1, 1, 1, 0, 1, 1] >>> [de() for i in range(10)] [4, 4, 5, 1, 6, 4, 3, 6, 2, 5]
Une expérience (une fonction, une réalisation), la répétition ($N$ appels, un compteur), la fréquence compteur/N. Au moins un $6$ en quatre lancers : $1-\left(\frac56\right)^4\simeq0{,}5177$.
def experience() : for i in range(4) : if de()==6 : return True return False def frequence(N) : compteur=0 for i in range(N) : if experience() : compteur=compteur+1 return compteur/N >>> random.seed(2026) >>> frequence(1000) 0.506 >>> frequence(100000) 0.5174 >>> 1-(5/6)**4 0.5177469135802468
$X\hookrightarrow\mathcal B(n,p)$ compte les succès de $n$ épreuves de Bernoulli indépendantes : on additionne $n$ appels à bernoulli(p), puis on compare les fréquences à $\mathbb P(X=k)=\binom nk p^k(1-p)^{n-k}$.
from math import comb def binomiale(n,p) : s=0 for i in range(n) : s=s+bernoulli(p) return s random.seed(2026) n,p,N=10,0.3,10000 E=[binomiale(n,p) for i in range(N)] for k in range(6) : freq=E.count(k)/N proba=comb(n,k)*p**k*(1-p)**(n-k) print(k, freq, round(proba,4)) # Sortie : k, fréquence, P(X=k) 0 0.0313 0.0282 1 0.122 0.1211 2 0.2294 0.2335 3 0.2594 0.2668 4 0.2002 0.2001 5 0.1093 0.1029
La fréquence des succès oscille d'abord, puis se resserre autour de $p=0{,}3$ sans s'y fixer : $\dfrac{S_n}{n}$ s'approche de $p$ avec une probabilité qui tend vers $1$.
def frequences(p,n) : F=[] s=0 for k in range(1,n+1) : s=s+bernoulli(p) F.append(s/k) return F >>> random.seed(2026) >>> F=frequences(0.3,2000) >>> F[9], F[99], F[999], F[1999] (0.3, 0.22, 0.286, 0.291)
from statistics import fmean, pstdev random.seed(2026) for n in [100,400,1600] : M=[binomiale(n,0.3)/n for i in range(1000)] print(n, round(fmean(M),4), round(pstdev(M),4), round(sqrt(0.3*0.7/n),4)) # Sortie : n, moyenne, écart-type observé, théorique 100 0.301 0.0478 0.0458 400 0.3012 0.0233 0.0229 1600 0.3006 0.0116 0.0115
$M_n=\dfrac{S_n}{n}$ a pour espérance $p$ et pour écart-type $\sigma(M_n)=\sqrt{\dfrac{p(1-p)}{n}}$ : multiplier $n$ par $4$ divise $\sigma$ par $2$, comme le confirme le tableau.
Pour $X\hookrightarrow\mathcal B(100\,;\,0{,}5)$, d'espérance $50$ et de variance $25$ :
random.seed(2026) n,p,N=100,0.5,10000 compteur=0 for i in range(N) : if abs(binomiale(n,p)-n*p)>=10 : compteur=compteur+1 print(compteur/N) print(2*sum(comb(100,k)/2**100 for k in range(41))) # Sortie : 0.0581 0.05688793364098079
La vraie valeur, environ $0{,}057$, est quatre fois plus petite : l'inégalité est grossière mais universelle.
Un point uniforme de $[0,1]^2$ tombe dans le quart de disque $x^2+y^2\leqslant1$ avec la probabilité $\dfrac\pi4$ : $4\times$ fréquence estime $\pi$.
def estim_pi(N) : compteur=0 for i in range(N) : x=random.random() y=random.random() if x**2+y**2<=1 : compteur=compteur+1 return 4*compteur/N >>> random.seed(2026) >>> estim_pi(10**6) 3.146604
L'erreur est de l'ordre de $\dfrac1{\sqrt N}$ (ici $\simeq\dfrac{1{,}64}{\sqrt N}$) : gagner une décimale exige de multiplier $N$ par $100$.
def aire_mc(f,a,b,M,N) : compteur=0 for i in range(N) : x=random.uniform(a,b) y=random.uniform(0,M) if y<=f(x) : compteur=compteur+1 return (b-a)*M*compteur/N def f(x) : return log(1+x) >>> random.seed(2026) >>> aire_mc(f,0,1,1,10**4) 0.389 >>> aire_mc(f,0,1,1,10**6) 0.385992
Un million de tirages donne trois décimales de $2\ln2-1\simeq0{,}386\,294$, un million de rectangles en donnait six : en dimension $1$, Monte-Carlo perd. Son intérêt : les domaines compliqués et la grande dimension.
def marche(n) : x=0 X=[0] for i in range(n) : x=x+random.choice([-1,1]) X.append(x) return X >>> random.seed(2026) >>> marche(15) [0, -1, 0, -1, -2, -1, 0, 1, 0, -1, -2, -3, -2, -3, -2, -3]
random.seed(2026) for n in [100,400,1600] : M=[marche(n)[-1] for i in range(2000)] print(n, fmean([abs(m) for m in M]), fmean([m*m for m in M])) # Sortie : n, moyenne de |X_n|, moyenne de X_n² 100 8.043 101.242 400 15.867 400.886 1600 32.929 1714.178
Pas de $\pm1$ équiprobables et indépendants : $\mathbb E(X_n)=0$ et $\operatorname{Var}(X_n)=n$. La distance typique à l'origine est de l'ordre de $\sqrt n$ : elle double quand $n$ est multiplié par $4$.
random()<p a pour probabilité $p$randint(1,6), bornes inclusesrandint(a,b) peut renvoyer $b$, contrairement à range(a,b).shuffle(L) modifie L et renvoie None : n'écrivez pas L=random.shuffle(L).seed.random.random()<p.