Mathématice, intégration des Tice dans l'enseignement des mathématiques  
Sommaire > Articles à paraître dans les numéros à venir > Mathématiques complémentaires : les algorithmes (...)

Mathématiques complémentaires : les algorithmes du programme de l’option de Terminale (2020).
Moteur de recherche
Mis en ligne le 23 mai 2020, par Benjamin Clerc

Cet article est long, complet et techniquement délicat : il rendra service à de très nombreux collègues.

Pourrions-nous demander aux utilisateurs de l’article de signaler les inévitables coquilles ou les anomalies qu’ils pourraient découvrir et qui nous auraient échappé...? Merci. (mathematice@sesamath.net)

Cet article peut être librement diffusé et son contenu réutilisé pour une utilisation non commerciale (contacter l’auteur pour une utilisation commerciale) suivant la licence CC-by-nc-sa

On peut lire dans le programme officiel de cette option :

Algorithmique et programmation
La démarche algorithmique est, depuis les origines, une composante essentielle de l’activité mathématique. Au collège, en mathématiques et en technologie, les élèves ont appris à écrire, mettre au point et exécuter un programme simple. Les classes de seconde et de première ont permis de consolider les acquis du collège (notion de variable, type de variables, affectation, instruction conditionnelle, boucle, notamment), d’introduire et d’utiliser la notion de fonction informatique et de liste. En algorithmique et programmation, le programme de mathématiques complémentaires reprend les programmes des classes de seconde et de première sans introduire de notion nouvelle, afin de consolider le travail des classes précédentes. Les algorithmes peuvent être écrits en langage naturel ou utiliser le langage Python. On utilise le symbole « ← » pour désigner l’affection dans un algorithme écrit en langage naturel. L’accent est mis sur la programmation modulaire qui permet de découper une tâche complexe en tâches plus simples. L’algorithmique trouve naturellement sa place dans toutes les parties du programme et aide à la compréhension et à la construction des notions mathématiques.

On va donc s’attacher à programmer les algorithmes listés dans le programme sans introduire de notions nouvelles :

Résolution d´équations du type $f(x) = k$ par balayage, par dichotomie, par la méthode de Newton.

Résolution d´équations par balayage
Cette méthode, classique, consiste, pour une fonction continue et monotone sur un intervalle $[a ; b]$ dans lequel on sait trouver la solution de $f(x) = k$, à calculer $f(x)$, pour $x$ allant de $a$ à $b$ avec un pas de $h$, après avoir vérifié la croissance ou la décroissance de la fonction. Si $f$ est croissante, "Tant Que" $f(x) < k$ est vérifiée (On utilisera pour cela un While.), on calcule $f(x)$. Dès que $f(x) > k$ on arrête l’algorithme. Si $f$est décroissante, "Tant Que" $f(x) > k$ est vérifiée, on calcule $f(x)$. Dès que $f(x) < k$ on arrête l’algorithme. On obtient un encadrement de longueur $h$ de la solution.
Il faut déterminer a et b, les bornes de l’intervalle dans lequel on va appliquer le balayage.
Fonction qui calcule cet encadrement, renvoyé dans un tuple, en fonction du nombre n de décimales souhaitées :
Cet algorithme détermine par balayage un encadrement de racine de 2 d’amplitude $10^{-n}$.

  1. from math import cos
  2. def f(x):
  3. return cos(x)-x
  4. def balayage(f,k,a,b,n):
  5. fa,fb=f(a),f(b)
  6. h=10**(-n)
  7. x=a
  8. if fa<fb:
  9. while f(x)<k:
  10. x=x+h
  11. else:
  12. while f(x)>k:
  13. x=x+h
  14. return (round(x-h,n),round(x,n))

Télécharger

Utilisation : Compiler le fichier puis dans la console :

  1. >> balayage(f,0,-1,2,6)
  2. (0.739085, 0.739086)

Télécharger

Le round permet ici de corriger quelques erreurs d’affichage des flottants [1].
Il est à noter que cet algorithme est très long à donner une réponse si n>=9.
On peut donc chercher à l’optimiser, par exemple en utilisant une boucle for qui va faire la même chose que le programme précédent mais pour h allant de 10-1 à 10-n en réduisant l’intervalle [a ; b] à chaque boucle.
C’est beaucoup plus rapide, jusqu’à n=16 où ça bogue, parce que les flottants (décimaux) ne sont affichés qu’avec 16 chiffres après la virgule.

  1. def balayage2(f,k,a,b,n):
  2. fa,fb=f(a),f(b)
  3. for i in range(1,n+1):
  4. h=10**(-i)
  5. x=a
  6. if fa<fb:
  7. while f(x)<k:
  8. x=x+h
  9. else:
  10. while f(x)>k:
  11. x=x+h
  12. a=x-h
  13. return (round(x-h,n),round(x,n))

Télécharger

Résolution d´équations par dichotomie

On va commencer par programmer l’algorithme classique de recherche de zéro d’une fonction sur un intervalle par dichotomie :

  1. def zero(f,a,b,epsilon):
  2. if f(a)*f(b)>0:
  3. return "Les bornes de l'intervalle ne sont pas bien choisies."
  4. while(abs(a-b)>epsilon):
  5. m=(a+b)/2
  6. if f(m)*f(a)>0:
  7. a=m
  8. else:
  9. b=m
  10. return m

Télécharger

Utilisation :

  1. >>> zero(f,-1,0,10**(-6))
  2. "Les bornes de l'intervalle ne sont pas bien choisies."
  3. >>> zero(f,-1,2,10**(-6))
  4. 0.7390849590301514

Télécharger

Et l’on peut aisément résoudre l’équation $f(x) = k$ avec un petit aménagement :

  1. def zero_k(f,a,b,epsilon,k):
  2. def fk(x):
  3. return f(x)-k
  4. return zero(fk,a,b,epsilon)

Télécharger

Utilisation :

  1. >>> zero_k(f,-20,0,10**(-4),10)
  2. -10.486984252929688

Télécharger

Résolution d´équations du type $f(x) = k$ par la méthode de Newton.

Soit $f$ une fonction dérivable sur un intervalle $I$.
L’équation $f (x) = 0$ admet une racine unique $\alpha$ sur cet intervalle $I$.
Soit $a \in I$ une valeur approchée de $\alpha$.
On va utiliser l’approximation affine $g$ de $f$ au point $a$.
On aura donc $g(x) = f (a) + (x − a)f ′(a)$ (tangente $T_a$).
La droite $T_a$ coupe l’axe des abscisses en $b = a −\frac{f(a)}{f’(a)}$.
Sous certaines conditions, le nombre $b$ peut représenter une meilleure approximation de $\alpha$ que $a$.
La méthode de Newton consiste à itérer le processus en repartant de $b$ et ainsi de suite.
On proposera deux solutions :
- si l’on sait dériver la fonction $f$.
  1. from math import sin
  2. def derivee_f(x):
  3. return -sin(x)-1
  4. def newton(x0,e):
  5. x = x0
  6. while abs(f(x))>e:
  7. d=derivee_f(x)
  8. x = x-f(x)/d
  9. return x

Télécharger


Utilisation :

  1. >>> newton(2,1e-6)
  2. 0.7390851332198145

Télécharger

PNG - 235.3 ko

Une lectrice avertie m’a fait une remarque qui m’amène à vous proposer une autre version :

  1. def newton2(fonction,derivee,x0,e):
  2. x = x0
  3. while abs(fonction(x))>e:
  4. d=derivee(x)
  5. x = x-fonction(x)/d
  6. return x

Télécharger

Utilisation :

  1. >>> newton2(f,derivee_f,2,1e-6)
  2. 0.7390851332198145

Télécharger


- si l’on doit se contenter d’un calcul approché de ce nombre dérivé en prenant comme approximation $f’(a) \approx \frac{f(a + h) - f(a)}{h}$ avec h suffisamment petit.

  1. def newton3(f,x0,e,h):
  2. x = x0
  3. while abs(f(x))>e:
  4. d = (f(x+h)-f(x))/h #on calcule une approximation du nombre dérivé en x
  5. x = x-f(x)/d
  6. return x

Télécharger

Utilisation :

  1. >>> newton3(f,2,1e-6,1e-5)
  2. 0.7390851332299205

Télécharger

Recherche d’une valeur approchée de précision donnée (Dans le thème d’étude Calculs d’aires)

Méthode des rectangles, des trapèzes.
On note $(x_i)_{i \in\{0, ..., n\}}$ une subdivision de l’intervalle [a ; b] Avec $a=x_0 < x_1 < ... < x_n = b$ et pour tout $i \in \{0, ..., n\}, x_i=a+i{\frac {b-a}{n}}$.
Il y a plusieurs méthodes pour approcher l’aire sous la courbe d’une fonction continue à l’aide de rectangles, commençons par la méthode des rectangles à gauche :
On remplace $f$ par la valeur qu’elle prend sur le bord gauche de l’intervalle $[x_i ; x_{i + 1}]$ soit $f(x_i)$, l’aire est alors égale à $f(x_i)\times{\frac{b - a}{n}}$.

D’où l’algorithme suivant :

  1. def f(x):
  2. return x**2-3*x+5
  3. def rectangles_gauche(f,a,b,n):
  4. h=(b-a)/n
  5. L=0
  6. for i in range(n):
  7. L+=f(a+i*h)
  8. return L*h

Télécharger

Utilisation :

  1. >>> rectangles_gauche(f,1,5,30)
  2. 24.54518518518518

Télécharger

Pour la méthode des rectangles à droite il suffit d’un tout petit changement :

  1. def rectangles_droit(f,a,b,n):
  2. h=(b-a)/n
  3. L=0
  4. for i in range(1,n+1):
  5. L+=f(a+i*h)
  6. return L*h

Télécharger

Utilisation :

  1. >>> rectangles_droit(f,1,5,30)
  2. 26.145185185185184

Télécharger

Ces méthodes minorent ou majorent le calcul d’aire selon les variations de $f$.
Une variante de ces deux méthodes consiste à prendre comme hauteur de chaque rectangle l’image du centre de l’intervalle d’où :

  1. def methode_milieux(f,a,b,n):
  2. h=(b-a)/n
  3. L=0
  4. for i in range(n):
  5. L+=f(a+(i+0.5)*h)
  6. return L*h

Télécharger

Utilisation :

  1. >>> methode_milieux(f,1,5,30)
  2. 25.327407407407406

Télécharger

Au lieu de remplacer $f$ par une fonction constante sur l’intervalle $[x_i ; x_{i + 1}]$, on peut la remplacer par une fonction affine représentée par le segment qui joint les points $(x_i ; f(x_{i}))$ et $(x_{i + 1} ; f(x_{i + 1}))$ :

  1. def trapezes(f,a,b,n):
  2. h=(b-a)/n
  3. L=0
  4. for i in range(n) :
  5. L+=(f(a+i*h)+f(a+i*h+h))/2
  6. return L*h

Télécharger

Utilisation :

  1. >>> trapezes(f,1,5,30)
  2. 25.345185185185183

Télécharger

Méthode de Monte-Carlo pour un calcul d’aire.

Cette méthode a été abordée en classe de première mais c’était limité au cas où $0\leq x\leq 1$ et $0\leq y\leq 1$, on va étendre la méthode au cas où $a\leq x\leq b$ et $0\leq y\leq M$ où M est un majorant de $f(x)$ sur $[a ; b]$.

On va tirer au hasard N points de coordonnées comprises entre a et b pour l’abscisse et entre 0 et M pour l’ordonnée, on compte le nombre de points S sous la courbe, l’aire sous la courbe est alors approchée par S/N*(b-a)*M.

  1. from random import *
  2. def f(x):
  3. return x**2-3*x+5
  4. def monte_carlo(f,a,b,M,n):
  5. S = 0
  6. for k in range(n):
  7. (x,y) = (uniform(a,b),uniform(0,M))
  8. if y < f(x):
  9. S += 1
  10. return S/n*M*(b-a)

Télécharger

Utilisation :

  1. >>> monte_carlo(f,1,5,25,1000000)
  2. 25.343799999999998

Télécharger

Chacune des méthodes précédentes peut être programmée en retournant directement la somme des éléments d’une liste plutôt que d’utiliser une boucle :

  1. def rectangles_gauche2(f,a,b,n):
  2. h=(b-a)/n
  3. return sum([f(a+i*h) for i in range(n)])*h
  4. def rectangles_droit2(f,a,b,n):
  5. h=(b-a)/n
  6. return sum([f(a+i*h) for i in range(1,n+1)])*h
  7. def methode_milieux2(f,a,b,n):
  8. h=(b-a)/n
  9. return sum([f(a+(i+0.5)*h) for i in range(n)])*h
  10. def trapezes2(f,a,b,n):
  11. h=(b-a)/n
  12. return sum([(f(a+i*h)+f(a+i*h+h))/2 for i in range(n)])*h
  13. def monte_carlo2(f, a, b, n):
  14. return sum([f(uniform(a, b)) for i in range(n)])*(b-a)/n

Télécharger

Pour cette dernière fonction, il n’a pas été possible de traduire la fonction monte_carlo() précédente à l’aide d’une liste, on a précédemment calculé une approximation de l’aire en procédant de la manière suivante :
- On a enfermé l’aire à calculer dans un rectangle de dimensions $(b - a)\times M$ ;
- On a lancé aléatoirement un grand nombre de points dans ce rectangle ;
- On a compté à chaque lancer ceux qui étaient sous la courbe, c’est-à-dire dans la surface que l’on veut calculer.
- On en a déduit par proportionnalité à une approximation de l’aire.
On ne peut pas simuler cela avec une liste mais on peut utiliser une variante de la méthode de Monte-Carlo (L’algorithme de la moyenne ou l’algorithme de l’espérance ou encore ici) :
- On tire aléatoirement et uniformément un grand nombre de nombres aléatoires x entre a et b ;
- On prend, à chaque tirage, comme approximation de l’aire, l’aire du rectangle de largeur b - a et de hauteur f(x) ;
- On calcule la moyenne des aires obtenues.
Il s’agit de l’algorithme monte_carlo2() ci-dessus.

Calcul des termes d’une suite - Pour une suite récurrente Un+1=f(Un), calcul des termes successifs - Calcul d’un terme de rang donné d’une suite

Les listes en compréhension sont des objets qui permettent de générer assez simplement des suites de nombres (voir Notion de liste en dernier paragraphe de cet article) :

  1. >>> [7*i for i in range(1,11)]
  2. [7, 14, 21, 28, 35, 42, 49, 56, 63, 70]
  3. >>> [2*i for i in range(20)]
  4. [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38]
  5. >>> [i**2 for i in range(16)]
  6. [0, 1, 4, 9, 16, 25, 36, 49, 64, 81, 100, 121, 144, 169, 196, 225]
  7. >>> [i for i in range(100) if i%7==2]
  8. [2, 9, 16, 23, 30, 37, 44, 51, 58, 65, 72, 79, 86, 93]

Télécharger

Voyons maintenant, plus classiquement, le calcul de termes d’une suite définie de manière explicite, c’est-à-dire à l’aide d’une fonction de n.

  1. def u(n):
  2. return 3*n+1
  3. def suite_explicite(u,debut,fin):
  4. suite = []
  5. for i in range(debut,fin+1):
  6. suite.append(u(i))
  7. return suite

Télécharger

Après avoir défini la fonction u, on crée une liste vide appelée "suite" dans laquelle on ajoute (append) à l’aide d’une boucle les termes d’indice compris entre début et fin.
Si l’on veut simplement un terme de rang donné, il suffit d’utiliser la fonction u :

  1. >>> suite_explicite(u,0,10)
  2. [1, 4, 7, 10, 13, 16, 19, 22, 25, 28, 31]
  3. >>> u(5)
  4. 16

Télécharger

Une autre version avec une liste :

  1. def suite_explicite2(u,debut,fin):
  2. return [u(i) for i in range(debut,fin+1)]

Télécharger

Regardons maintenant du côté des suites définies par récurrence :

  1. def f(u_n):
  2. return 3*u_n+1
  3. def suite_recurrente(f,debut,fin,u_0):
  4. suite = [u_0]
  5. for i in range(debut,fin):
  6. suite.append(f(suite[i-debut]))
  7. return suite

Télécharger

Après avoir défini la fonction f, qui calcule $u_{n+1}$ en fonction de $u_{n}$, on crée une liste appelée "suite" qui contient uniquement le premier terme de la liste, dans laquelle on ajoute (append) à l’aide d’une boucle les termes d’indice compris entre début+1 et fin.

  1. >>> suite_recurrente(f,0,10,0)
  2. [0, 1, 4, 13, 40, 121, 364, 1093, 3280, 9841, 29524]
  3. >>> suite_recurrente(f,5,10,12)
  4. [12, 37, 112, 337, 1012, 3037]

Télécharger

Utilisation :

  1. >>> suite_recurrente(f,0,10,0)
  2. 29524
  3. >>> suite_recurrente(f,5,10,12)
  4. 3037

Télécharger

Si l’on veut simplement un terme de rang "fin", il suffit de renvoyer suite[fin-debut] à la place de suite.

  1. def suite_recurrente(f,debut,fin,u_0):
  2. suite = [u_0]
  3. for i in range(debut,fin):
  4. suite.append(f(suite[i-debut]))
  5. return suite[fin-debut]

Télécharger

Utilisation :

  1. >>> suite_recurrente(f,5,10,12)
  2. 3037
  3. >>> suite_recurrente(f,0,10,0)
  4. 29524

Télécharger

Recherche de seuils

Un algorithme "de seuil" permet de déterminer une valeur pour laquelle une condition est respectée pour la première fois.
Ce type d’algorithme peut être utilisé pour répondre, par exemple, à ce genre de question : « déterminer la plus petite valeur de n pour laquelle $u_n>1000$ ».

Concrètement, on considère une suite $(u_n)$ et on cherche à écrire un algorithme "de seuil" sur cette suite.
Le modèle à utiliser est le suivant (où les termes de la suite sont représentés par la lettre U et les rangs par la lettre N) :

N ← valeur initiale (souvent 0 ou 1)
U  ← valeur initiale
TANT QUE condition
   N  ← N+1
   U  ← .... (le terme suivant qui dépend de l'exercice)
FIN TANT QUE

On en profite pour calculer la somme des termes de la suite entre les indices debut et fin.
Si la suite est définie de manière explicite :

  1. def u(n):
  2. return 3*n+1
  3. def seuil_explicite(n_0,N):
  4. n=n_0
  5. while u(n)<N:
  6. n=n+1
  7. return n
  8. def somme_explicite(debut,fin):
  9. s=0
  10. for i in range(debut,fin+1):
  11. s=s+u(i)
  12. return s

Télécharger

Utilisation :

  1. >>> seuil_explicite(5,1000)
  2. 333
  3. >>> somme_explicite(5,9)
  4. 110

Télécharger

Si la suite est définie de manière récurrente :

  1. def f(u_n):
  2. return 3*u_n+1
  3. def seuil_recurrent(n_0,u_0,N):
  4. n=n_0
  5. U=u_0
  6. while U<N:
  7. n=n+1
  8. U=f(U)
  9. return n
  10. def somme_recurrente(debut,fin,u_0):
  11. s=u_0
  12. u=u_0
  13. for i in range(debut,fin):
  14. u=f(u)
  15. s=s+u
  16. return s

Télécharger

Utilisation :

  1. >>> seuil_recurrent(5,1000,15000)
  2. 8
  3. >>> somme_recurrente(3,7,5)
  4. 663

Télécharger

On pourra facilement adapter ces scripts de seuil si la suite est décroissante et qu’il s’agit de déterminer à partir de quand son terme devient inférieur à une valeur donnée.

Recherche de valeurs approchées de constantes mathématiques, par exemple π, ln2,√2

Valeur approchée de π
Nous avons déjà approximé π en classe de première avec la méthode d’Archimède et la méthode de Monte Carlo, nous utiliserons en Terminale la formule de Bellard :
$\displaystyle \pi ={\frac {1}{2^{6}}}\sum _{n=0}^{\infty }{\frac {(-1)^{n}}{2^{10n}}}\left(-{\frac {2^{5}}{4n+1}}-{\frac {1}{4n+3}}+{\frac {2^{8}}{10n+1}}-{\frac {2^{6}}{10n+3}}-{\frac {2^{2}}{10n+5}}-{\frac {2^{2}}{10n+7}}+{\frac {1}{10n+9}}\right)$

  1. def bellard(n):
  2. pi=0
  3. for i in range(n+1):
  4. pi+=(-2**5/(4*i+1)-1/(4*i+3)+2**8/(10*i+1)-2**6/(10*i+3)-2**2/(10*i+5)-2**2/(10*i+7)+1/(10*i+9)) *(-1)**i/(2**6*2**(10*i))
  5. return pi

Télécharger

Utilisation :

  1. >>> bellard(3)
  2. 3.1415926535897545
  3. >>> bellard(4)
  4. 3.1415926535897922

Télécharger

A noter que Python ne peut faire mieux que bellard(4) qui n’est pas une bonne approximation puisque les 3 dernières décimales sont erronées. En utilisant le module Decimal, on obtient deux décimales exactes de plus, pas plus ...

Valeur approchée de $\ln{2}$
En écrivant $\ln(2)=\int_{1}^{2} \frac{1}{t}dt$ et en subdivisant de manière de plus en plus fine l’intervalle [1 ; 2] pour approcher cette intégrale par une somme de Riemann.

Nous verrons cela plus loin dans l’article.
Nous allons donc calculer $ln(2)$ au moyen de la série harmonique alternée : ${\displaystyle \ln(2) =\lim _{n\to \infty }\left(\sum _{k=1}^{n}{\frac {(-1)^{k + 1}}{k}}\right)}$

  1. def ln2(n):
  2. ln2=0
  3. for i in range(1,n+1):
  4. ln2+=(-1)**(i+1)/i
  5. return ln2

Télécharger

Utilisation :

  1. >>> ln2(1000)
  2. 0.6926474305598223
  3. >>> ln2(10000)
  4. 0.6930971830599583

Télécharger

Cela converge très lentement puisque $\ln(2) \approx 0.69314718056$ ...

Valeur approchée de $\sqrt{2}$
On doit à Théon de Smyrne ces deux suites $(p_n)$ et $(q_n)$ définies par récurrence :
$p_{n + 1} = p_n + 2q_n, p_0 = 1$ ;
$q_{n + 1} = p_n + q_n, q_0 = 1$.
Ces suites sont à valeur entière strictement positive, donc strictement croissantes par récurrence, et vérifient :
$p_n² − 2q_n² = (−1)^n(p_0² − 2q_0²)$ de sorte que $\frac{p_n}{q_n}$ tend vers $\sqrt{2}$.

  1. def racine2(n):
  2. p,q=1,1
  3. for i in range(n):
  4. p,q=p+2*q,p+q
  5. return p/q

Télécharger

Utilisation :

  1. >>> racine2(10)
  2. 1.4142135516460548
  3. >>> racine2(100)
  4. 1.4142135623730951

Télécharger

Méthode d’Euler

Soit $f$ une fonction dérivable sur $]0 ; +\infty[$ et vérifiant $f(1) = 0$ et pour tout $x \in ]0 ; +\infty[, f’(x) = \frac{1}{x}$. Utilisons la méthode d’Euler pour tracer une courbe approchée sur l’intervalle [1 ; 3] de cette fonction f.
La méthode d’Euler consiste à construire une suite de points $(x_n ; y_n)$ telle que $y_n$ soit proche de $f(x_n)$. Ainsi le nuage de points $(x_n ; y_n)$ formera une approximation de la courbe représentative de la fonction $f$.
La suite $(x_n)$ est une suite arithmétique de raison $h>0$ et de premier terme 1, $h$ est appelé le pas. On a donc $x_n = 1 + nh$ et $x_{n+1}=x_n + h$.
La suite $(y_n)$ est définie par $y_0 = f(x_0) = f(1) = 0$ et pour tout $n \in \mathbf{N}, y_{n+1}=y_n+hf’(x_n)$. C’est-à-dire que l’on se sert de l’approximation affine de la fonction par ses tangentes successives.
On a $y_{n+1}=y_n+\frac{h}{x_n}$.
D’où le script suivant :

  1. import matplotlib.pyplot as plt
  2. def Euler(x_f,n) :
  3. x=1
  4. y=0
  5. h=x_f/n
  6. abscisse=[1]
  7. ordonnee=[0]
  8. for i in range(n) :
  9. y=y+h/x
  10. x=x+h
  11. abscisse.append(x)
  12. ordonnee.append(y)
  13. plt.plot(abscisse,ordonnee,"r")
  14. plt.show()

Télécharger

Utilisation, pour un nuage de n=100 points pour x allant de 1 à x_f=3 :

  1. >>> Euler(3,100)

Ce qui donne le graphique suivant :

Et si l’on veut comparer avec le "vrai" graphe de la fonction logarithme neperien :

  1. import matplotlib.pyplot as plt
  2. from math import log
  3. def Euler(x_f,n) :
  4. x=1
  5. y=0
  6. h=x_f/n
  7. abscisse=[1]
  8. ordonnee=[0]
  9. ln=[0]
  10. for i in range(n) :
  11. y=y+h/x
  12. x=x+h
  13. abscisse.append(x)
  14. ordonnee.append(y)
  15. ln.append(ln(x))
  16. plt.plot(abscisse,ordonnee,"r")
  17. plt.plot(abscisse,ln,"b")
  18. plt.show()

Télécharger

Algorithme de Briggs pour le calcul de logarithmes

L’Algorithme de Briggs est basé sur le fait que $\ln(\sqrt{x}) = \frac 1 2 \ln(x)$ (Voir page 50 de ce document du site exo7 et l’article d’Alain Busser).

Entrée : Un réel x strictement positif
Variable : Un entier n
Tant que x est éloigné de 1 (à 0,001 près par exemple), faire
       incrémenter le compteur n ;
       remplacer x par sa racine carrée
#Comme le logarithme népérien d’un nombre x proche de 1 est proche de x-1,
#on prend x-1 comme valeur approchée du logarithme ; n vaut alors le nombre de fois qu’on doit doubler x-1 pour avoir le logarithme du nombre de départ.
Faire n fois :
       Multiplier le logarithme par 2
Sortie : Le nombre obtenu est le logarithme népérien de x.

Ce qui donne en Python, en prenant un second paramètre p pour arrêter la boucle Tant que dès que x est éloigné de 1 de moins de $10^{-p}$ :

  1. from math import *
  2. def ln(x,p):
  3. n=0
  4. while(abs(x-1)>10**(-p)):
  5. n+=1
  6. x=sqrt(x)
  7. y=x-1
  8. while(n>0):
  9. y=y*2
  10. n=n-1
  11. return y

Télécharger

Utilisation :

  1. >>> ln(5,3)
  2. 1.6100704732402846
  3. >>> ln(5,8)
  4. 1.609437882900238
  5. >>> from math import log
  6. >>> log(5)
  7. 1.6094379124341003

Télécharger

Alors que la calculatrice donne 1.609437912 et WolframAlpha 1.609437912434100374600759333226187639525601354268517721912...
On constate que la précision pour ln(5,8) est donc d’environ $3\times 10^{-8}$

Approximation de ln 2 par dichotomie selon l’algorithme de Brouncker

On a précédemment approché $ln(2)$ avec la formule $ln(2) \approx 1 - \frac{1}{2} + \frac{1}{3} - \frac{1}{4} + \frac{1}{5} - \frac{1}{6} + \frac{1}{7} ...$ or on peut remarquer que $1 - \frac{1}{2} = \frac{1}{1\times 2}$, $\frac{1}{3} - \frac{1}{4} = \frac{1}{3\times 4}$, $\frac{1}{5} - \frac{1}{6} = \frac{1}{5\times 6}$, ... donc on peut écrire $ln(2)\approx \frac{1}{1\times 2} + \frac{1}{3\times 4} + \frac{1}{5\times 6} + + \frac{1}{7\times 8} + ...$ d’où la formule de Lord Brouncker ${\displaystyle \ln(2) =\lim_{n\to \infty} \left(\sum_{k=0}^{n}{\frac {1}{(2k+1)(2k+2)}}\right)}$ d’où l’algorithme suivant :

  1. def brouncker(n):
  2. ln2=0
  3. for i in range(n):
  4. ln2+=1/((2*i+1)*(2*i+2))
  5. return ln2

Télécharger

Utilisation :

  1. >>> brouncker(20)
  2. 0.6808033817926938
  3. >>> brouncker(200)
  4. 0.691898743055063
  5. >>> from math import log
  6. >>> log(2)
  7. 0.6931471805599453

Télécharger

On peut voir que ça ne converge pas très vite même si cela est plus rapide qu’avec la série alternée.

Dans le cadre de la loi binomiale : calcul de coefficients binomiaux (triangle de Pascal), de probabilités ; détermination d’un intervalle I pour lequel la probabilité P(XI) est inférieure à une valeur donnée $\alpha$, ou supérieure à $1 - \alpha$.

Calcul de coefficients binomiaux (triangle de Pascal)

Le triangle de Pascal :

1
1 1
1 2 1
1 3 3 1
1 4 6 4 1
1 5 10 10 5 1
$\begin{pmatrix}0\\0\end{pmatrix}$
$\begin{pmatrix}1\\0\end{pmatrix} \begin{pmatrix}1\\1\end{pmatrix}$
$\begin{pmatrix}2\\0\end{pmatrix}\begin{pmatrix}2\\1\end{pmatrix}\begin{pmatrix}2\\2\end{pmatrix}$
$\begin{pmatrix}3\\0\end{pmatrix}\begin{pmatrix}3\\1\end{pmatrix}\begin{pmatrix}3\\2\end{pmatrix}\begin{pmatrix}3\\3\end{pmatrix}$
$\begin{pmatrix}4\\0\end{pmatrix}\begin{pmatrix}4\\1\end{pmatrix}\begin{pmatrix}4\\2\end{pmatrix}\begin{pmatrix}4\\3\end{pmatrix}\begin{pmatrix}4\\4\end{pmatrix}$
$\begin{pmatrix}5\\0\end{pmatrix}\begin{pmatrix}5\\1\end{pmatrix}\begin{pmatrix}5\\2\end{pmatrix}\begin{pmatrix}5\\3\end{pmatrix}\begin{pmatrix}5\\4\end{pmatrix}\begin{pmatrix}5\\5\end{pmatrix}$

La construction du triangle suit la relation de Pascal :
$\begin{pmatrix}n\\k\end{pmatrix} = \begin{pmatrix}n - 1 \\k - 1 \end{pmatrix} + \begin{pmatrix}n - 1\\k\end{pmatrix}$
Python possède une commande binomial du module Sympy (Installation du module Sympy en installant Anaconda) qui permet de calculer le coefficient binomial de deux entiers naturels :

  1. >>> from sympy import binomial
  2. >>> binomial(6,2)
  3. 15

Télécharger

La commande comb() class="python" a récemment intégré le module math, on peut donc également utiliser :

  1. >>> from math import comb
  2. >>> comb(10,2)
  3. 45

Télécharger

Nous allons bien sûr nous passer de cela et programmer une fonction qui calcule $\begin{pmatrix}n\\k\end{pmatrix}$ par itération.
Au lieu de construire le triangle de Pascal en entier, nous pouvons nous focaliser sur la construction de la k-ième colonne en constatant que :
$\begin{pmatrix}n\\k + 1\end{pmatrix} =\frac{{n} !}{{(k + 1) !}{(n - (k + 1)) !}} = \frac{n !}{k ! (n - k) !}\frac{n - k}{k + 1} = \frac{n - k}{k + 1}\begin{pmatrix}n\\k\end{pmatrix} $
On peut donc générer les termes en itérant à partir de $\begin{pmatrix}1\\0\end{pmatrix}$

  1. def pascal(n,k):
  2. coeff = 1
  3. for i in range(k):
  4. coeff = (coeff*(n-i))//(i+1)
  5. return coeff

Télécharger

Utilisation :

  1. >>> pascal(6,2)
  2. 15
  3. >>> pascal(10,2)
  4. 45

Télécharger

Enfin si on utilise la commande factorial() du module math :

  1. from math import factorial
  2. def C(n,k):
  3. return factorial(n)//(factorial(k)*factorial(n-k))

Télécharger

Utilisation :

  1. >>> C(10,2)
  2. 45
  3. >>> C(6,2)
  4. 15

Télécharger

Simulation avec Python d’une variable aléatoire (de la loi de Bernoulli, d’une loi uniforme discrète, etc.) d’un échantillon de taille 𝑛 d’une variable aléatoire.

Loi de Bernoulli
La distribution de Bernoulli ou loi de Bernoulli, du nom du mathématicien suisse Jacques Bernoulli, est une distribution discrète de probabilité, qui prend la valeur 1 avec la probabilité p et 0 avec la probabilité q = 1 – p. D’où la fonction suivante :

  1. from random import random
  2. def epreuve(p):
  3. if random()<=p:
  4. return True
  5. else:
  6. return False

Télécharger

Un schéma de Bernoulli consiste à effectuer un certain de nombre de fois une même épreuve de Bernoulli, le résultat de chaque épreuve de Bernoulli étant indépendant des résultats des autres épreuves de Bernoulli. Le nombre d’épreuves est traditionnellement noté n :

  1. def schema(n,p):
  2. nb_succes=0
  3. for i in range(n):
  4. if epreuve(p)==True:
  5. nb_succes+=1
  6. return nb_succes

Télécharger

On considère alors X la variable aléatoire qui associe à une réalisation d’un schéma de Bernoulli le nombre de succès en n tentatives. On dit que la variable aléatoire X est régie par une loi binomiale de paramètres n et p. On note B(n,p) la loi binomiale de paramètres n et p. On peut alors réaliser une expérience qui consiste à répéter N fois le schéma de Bernoulli :

  1. def experience(N,n,p):
  2. frequences=[]
  3. for i in range(N):
  4. frequences.append(0)
  5. for i in range(N):
  6. nb_succes=schema(n,p)
  7. frequences[i]=nb_succes/n
  8. return frequences

Télécharger

Utilisation :

  1. >>> experience(10,1000,0.36)
  2. [0.354, 0.346, 0.329, 0.329, 0.39, 0.367, 0.368, 0.341, 0.352, 0.327]

Télécharger

Loi uniforme discrète
Pour une variable aléatoire X qui suit une loi discrète uniforme sur $\{1,2,...,k\}$, on va partitionner l’intervalle [0,1] dans lequel une variable aléatoire U qui suit une loi uniforme prend ses valeurs en k sous-intervalles de taille $\frac 1 k$. Il faut ensuite déterminer dans lequel de ces sous-intervalles se trouve U. Une manière de faire cela est d’utiliser la partie entière supérieure :

  1. from random import random
  2. from math import ceil
  3. def uniforme_discrete(k):
  4. U=random()
  5. X=ceil(k*U)
  6. return X

Télécharger

Utilisation :

  1. >>> for i in range(100):
  2. uniforme_discrete(10)
  3. 9
  4. 3
  5. 3
  6. 1
  7. 10
  8. 7
  9. 1
  10. 5
  11. 8
  12. 6

Télécharger

Et si l’on veut un échantillon de taille n :

  1. def schema(n,k):
  2. L=[]
  3. for i in range(k):
  4. L.append(0)
  5. for i in range(n):
  6. sortie=uniforme_discrete(k)
  7. for j in range(k):
  8. if sortie==j+1:
  9. L[j]+=1
  10. return L

Télécharger

Utilisation :

  1. >>> schema(100,10)
  2. [14, 14, 11, 5, 10, 8, 13, 9, 10, 6]
  3. >>> schema(1000,12)
  4. [81, 74, 80, 92, 82, 74, 91, 91, 78, 96, 77, 84]

Télécharger

Et si l’on veut réaliser l’expérience N fois :

  1. def experience(N,n,k):
  2. frequences=[]
  3. for i in range(k):
  4. frequences.append(0)
  5. for i in range(k):
  6. nb_succes=schema(n,k)
  7. frequences[i]=nb_succes[i]/n
  8. return frequences

Télécharger

Utilisation :

  1. >>> experience(10,1000,8)
  2. [0.12, 0.124, 0.124, 0.153, 0.136, 0.141, 0.13, 0.124]
  3. >>> experience(10,100000,8)
  4. [0.12371, 0.12461, 0.12588, 0.12389, 0.12671, 0.1243, 0.12541, 0.12683]

Télécharger

Loi discrète
Par exemple : Loi sur $\{-10,3,50\}$ de loi $P(-10) =\frac 2 3, P(3) =\frac 1 6, P(50) =\frac 1 6$. On partitionne l’intervalle [0,1] en trois sous-intervalles de taille respective $\frac 2 3, \frac 1 6 et \frac 1 6$. Une manière de faire consiste à prendre $[0 ; \frac 2 3], [\frac 2 3 ; \frac 2 3 + \frac 1 6], [\frac 2 3 + \frac 1 6 ; 1]$. L’algorithme est alors :

  1. from random import random
  2. def loi_discrete():
  3. U=random()
  4. if U<2/3:
  5. X=-10
  6. elif U<5/6:
  7. X= 3
  8. else:
  9. X= 50
  10. return X

Télécharger

Et pour simuler n expériences :

  1. def simule(n):
  2. L=[0,0,0]
  3. for i in range(n):
  4. sortie=loi_discrete()
  5. if sortie==-10:
  6. L[0]+=1
  7. elif sortie==3:
  8. L[1]+=1
  9. else:
  10. L[2]+=1
  11. return L

Télécharger

Utilisation :

  1. >>> simule(10000)
  2. [6650, 1708, 1642]

Télécharger

Et si l’on préfère récupérer les fréquences :

  1. for i in range(3):
  2. L[i]/=n

Télécharger

  1. >>> simule(10000)
  2. [0.6628, 0.1644, 0.1728]

Télécharger

Avec la commande choices() du module random :
Avec les poids cumulatifs :
Pour simuler la loi discrète :

  1. >>> from random import *
  2. >>> choices([-10,3,50],cum_weights=[2/3,5/6,1],k=1)
  3. [3]

Télécharger

Et pour simuler k expériences :

  1. >>> choices([-10,3,50],cum_weights=[2/3,5/6,1],k=15)
  2. [-10, 50, -10, 50, 50, -10, 3, -10, -10, -10, -10, 50, -10, 3, -10]

Télécharger

Si l’on veut récupérer le tableau des effectifs :

  1. def loi_discrete_choices(X,P,n):
  2. L=choices(X,cum_weights=P,k=n)
  3. l=len(X)
  4. S=[0]*l
  5. for i in range(l):
  6. S[i]=L.count(X[i])
  7. return S

Télécharger

Utilisation :

  1. >>> loi_discrete_choices([-10,3,50],[2/3,5/6,1],1000)
  2. [663, 180, 157]

Télécharger

Et si l’on veut les fréquences, on remplace S[i]=L.count(X[i]) par S[i]=L.count(X[i])/n
Utilisation :

  1. >>> loi_discrete_choices([-10,3,50],[2/3,5/6,1],1000)
  2. [0.694, 0.15, 0.156]

Télécharger

Avec les poids initiaux :
Pour simuler la loi discrète :

  1. >>> from random import *
  2. >>> choices([-10,3,50],[2/3,1/6,1/6],k=1)
  3. [3]

Télécharger

Et pour simuler n expériences :

  1. >>> choices([-10,3,50],[2/3,1/6,1/6],k=15)
  2. [3, 3, -10, 50, 3, -10, 3, 3, -10, 3, 50, -10, -10, 50, -10]

Télécharger

Si l’on veut récupérer le tableau des effectifs :

  1. def loi_discrete_choices(X,P,n):
  2. L=choices(X,P,k=n)
  3. l=len(X)
  4. S=[0]*l
  5. for i in range(l):
  6. S[i]=L.count(X[i])
  7. return S

Télécharger

Utilisation :

  1. >>> loi_discrete_choices([-10,3,50],[2/3,1/6,1/6],1000)
  2. [663, 180, 157]

Télécharger

Et si l’on veut les fréquences, on remplace S[i]=L.count(X[i]) par S[i]=L.count(X[i])/n
Utilisation :

  1. >>> loi_discrete_choices([-10,3,50],[2/3,1/6,1/6],1000)
  2. [0.694, 0.15, 0.156]

Télécharger

Fonction Python renvoyant une moyenne pour un échantillon. Série des moyennes pour N échantillons de taille n d’une variable aléatoire d’espérance μ et d’écart type σ. Calcul de l’écart type s de la série des moyennes des échantillons observés,à comparer à σ/√n. Calcul de la proportion des cas où l’écart entre m et μ est inférieur ou égal à kσ/√n ou à ks, pour k=2 ou k=3.

On simule des échantillons de taille N d’une variable aléatoire X suivant une loi discrète uniforme dans $\{1 ; 2 ; ... ; n\}$ (cela revient à lancer un dé équilibré à n faces) :

  1. from random import *
  2. def echantillon(N,n):
  3. L=[]
  4. for i in range(N):
  5. L.append(randint(1,n))
  6. return L

Télécharger

Utilisation :

  1. >>> echantillon(10,8)
  2. [6, 4, 2, 2, 4, 1, 8, 8, 7, 8]

Télécharger

On peut alors s’intéresser à la moyenne de cet échantillon :

  1. def moyenne(L):
  2. return sum(L)/len(L)

Télécharger

Et constater qu’en général la moyenne obtenue n’est pas égale à l’espérance :

  1. >>> moyenne(echantillon(10,8))
  2. 3.9

Télécharger

$E(X)=\frac{n+1}{2}$.
Dans notre exemple $E(X) = \frac{9}{2} = 4,5 \neq 3,9$.
Pour comprendre ce phénomène, on peut simuler des échantillons de taille n et calculer la moyenne de ces échantillons :

  1. def SerieMoyennes(NS,N,n):
  2. M=[]
  3. for i in range(NS):
  4. M.append(moyenne(echantillon(N,n)))
  5. return M

Télécharger

Utilisation :

  1. >>> SerieMoyennes(20,1000,8)
  2. [4.468, 4.526, 4.365, 4.605, 4.512, 4.536, 4.538, 4.534, 4.547, 4.462, 4.418, 4.493, 4.534, 4.509, 4.466, 4.493, 4.523, 4.457, 4.324, 4.338]

Télécharger

On n’obtient pas une seule fois 4,5 sur ces 20 séries mais la taille des échantillons étant assez grande, on trouve des moyennes "proches" de 4,5.
On va alors s’intéresser à cette notion de proximité en calculant l’écart-type :

  1. from math import *
  2. def EcartType(L):
  3. n=len(L)
  4. m=moyenne(L)
  5. ET=[]
  6. for i in range(n):
  7. ET.append((L[i]-m)**2)
  8. v=moyenne(ET)
  9. s=sqrt(v)
  10. return s

Télécharger

Puis la proportion de ces moyennes m qui appartiennent aux intervalles $[\mu - \frac{k\sigma}{\sqrt{N}} ; \mu + \frac{k\sigma}{\sqrt{N}}]$ :

  1. def prop_int_ksigma(NS,N,n,k):
  2. e=0
  3. S=SerieMoyennes(NS,N,n)
  4. mu=(1+n)/2
  5. sigma=EcartType(S)
  6. for m in S:
  7. if abs(m-mu)<k*sigma/sqrt(N):
  8. e=e+1
  9. return e/N

Télécharger

Utilisation :

  1. >>> prop_int_ksigma(1000,500,8,1)
  2. 0.084
  3. >>> prop_int_ksigma(1000,500,8,2)
  4. 0.118
  5. >>> prop_int_ksigma(1000,500,8,3)
  6. 0.192

Télécharger

Simulation d’une variable aléatoire de loi géométrique à partir du schéma de Bernoulli.

Définition : Soit $0 < p < 1$. On dit qu’une variable $X$ suit une loi géométrique de paramètre $p$, notée $G(p)$ si
a) $X$ est à valeurs dans $\mathbb{N}^{{}∗{}}$ ;
b) $∀k≥1, P(N=k) = (1−p)^{k−1}p$.
La variable aléatoire qui suit une loi géométrique de paramètre $p$ peut être assimilée au temps d’attente du premier succès d’un schéma de Bernoulli :

  1. from random import random
  2. def attente_succes(p):
  3. temps_attente=1
  4. alea=random()
  5. while alea>p:
  6. alea=random()
  7. temps_attente+=1
  8. return temps_attente

Télécharger

Utilisation :

  1. >>> attente_succes(0.2)
  2. 3
  3. >>> attente_succes(0.2)
  4. 10
  5. >>> attente_succes(0.2)
  6. 1
  7. >>> attente_succes(0.2)
  8. 4
  9. >>> attente_succes(0.2)
  10. 23

Télécharger

On peut alors simuler un grand nombre de fois l’expérience et représenter graphiquement le résultat :

  1. from random import random
  2. import matplotlib.pyplot as pyplot
  3. def attente_succes(p):
  4. temps_attente=1
  5. alea=random()
  6. while alea>p:
  7. alea=random()
  8. temps_attente+=1
  9. return temps_attente
  10. def simulation(p,N):
  11. liste_nbtirages=[0]*30
  12. for i in range(N):
  13. t=attente_succes(p)
  14. if t<=30:
  15. liste_nbtirages[t-1]+=1
  16. pyplot.xlabel("Nombde de tentatives avant le succès")
  17. pyplot.ylabel("Effectifs")
  18. pyplot.bar(range(1,31),liste_nbtirages,width=0.5,color="blue")
  19. pyplot.show()

Télécharger

Utilisation :

  1. >>> simulation(0.2,100000)


On peut aussi afficher le résultat théorique :

  1. def loi_geom(p,k):
  2. K=[k for k in range(1,k+1)]
  3. P=[p*(1-p)**(k-1) for k in range(1,k+1)]
  4. pyplot.xlabel("Nombre de tentatives avant le succès")
  5. pyplot.ylabel("Probabilité")
  6. pyplot.bar(K,P)
  7. pyplot.show()

Télécharger

Utilisation :

  1. >>> loi_geom(0.2,30)

Simulation d’une loi exponentielle à partir d’une loi uniforme.

On peut simuler la loi uniforme sur l’intervalle [a ; b] à l’aide de la fonction random() du module random :

  1. from random import random
  2. def uniforme(a,b):
  3. return(a+(b-a)*random())

Télécharger

En effet random() renvoie un nombre aléatoire de l’intervalle [0 ; 1[ (La suite de nombres générée “ressemble” à une suite de variables aléatoires uniformes et indépendantes.), donc (b-a)*random() renvoie un nombre aléatoire de l’intervalle [0 ; b-a[ et a+(b-a)*random() renvoie un nombre aléatoire de l’intervalle [a ; b[ .
La fonction densité de probabilité de cette loi uniforme sur [a ; b] est une fonction constante par morceaux :

  1. def densuni(x,a,b):
  2. #densité sur [a,b]
  3. if (a<=x and x<=b):
  4. return(1/(b-a))
  5. else:
  6. return(0)

Télécharger

Et si l’on veut sa représentation graphique :

  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3. from random import random
  4. def uniforme(a,b):
  5. return(a+(b-a)*random())
  6. def densuni(x,a,b):
  7. #densité sur [a,b]
  8. if (a<=x and x<=b):
  9. return(1/(b-a))
  10. else:
  11. return(0)
  12. def graph_uniforme(a,b):
  13. abs=np.linspace(a-(b-a),b+(b-a),200)
  14. y=[densuni(x,a,b) for x in abs]
  15. plt.grid()
  16. plt.plot(abs,y)
  17. plt.show()

Télécharger

Utilisation :

  1. >>>graph_uniforme(-5,10)


On peut aussi s’intéresser à la fonction de répartition :

  1. def repartuni(x,a,b):
  2. #fonction de répartition sur [a,b]
  3. if x<a:
  4. return(0)
  5. elif (a<=x and x<=b):
  6. return((x-a)/(b-a))
  7. else:
  8. return(1)
  9. def graph_repartuniforme(a,b):
  10. abs=np.linspace(a-(b-a),b+(b-a),200)
  11. y=[repartuni(x,a,b) for x in abs]
  12. plt.grid()
  13. plt.plot(abs,y)
  14. plt.show()

Télécharger

Utilisation :

  1. >>>graph_repartuniforme(-5,10)

Pour la loi exponentielle de paramètre λ, on a $F(x) = 1−e^{−λx}$ pour $x \geqslant 0$ et $F(x) = 0$ si $x \lt 0$. Sur $\mathbb{R}+$ , la fonction F(x) est continue et strictement croissante, on a l’équivalence :
$u = 1−e^{−λx} \Leftrightarrow 1 - u = e^{−λx} \Leftrightarrow ln(1 - u) = −λx \Leftrightarrow x = F^{−1}(u) =−\frac{1 }{λ}ln(1−u)$.
On peut donc simuler X à partir de U puisque $F^{−1}(U)$ suit la même loi que X.

On peut alors écrire une fonction Expo qui prend en argument un réel a>0 et renvoie une réalisation d’une variable aléatoire qui suit la loi Exponentielle de paramètre a :

  1. from math import log
  2. from random import *
  3. def Expo(a):
  4. return -1/a*log(1-random())

Télécharger

Utilisation :

  1. >>>Expo(2)
  2. >>>0.5242804058215232
  3. >>>Expo(0.2)
  4. >>>3.2637966061217183

Télécharger

On peut alors utiliser une fonction qui va procéder à un nombre n de tirages aléatoires :

  1. def tirages(n,a):
  2. L=[]
  3. for i in range(n):
  4. L.append(Expo(a))
  5. return L

Télécharger

Utilisation :

  1. >>>tirages(10,0.3)
  2. >>>[1.4373901275636258,
  3. 4.961359207330598,
  4. 1.5097403973808872,
  5. 5.805058647163496,
  6. 9.818932947901471,
  7. 11.037823931384885,
  8. 1.3647856789399562,
  9. 3.4588250161094725,
  10. 2.9218462025078207,
  11. 4.015415133657791]

Télécharger

On peut alors passer à la représentation graphique en utilisant la commande hist [2] du module matplotlib :

  1. def graphe(n,a):
  2. L=tirages(n,a)
  3. maxL=int(max(L))
  4. hist(L , bins = 100, range=(0,maxL), normed = True)

Télécharger


On peut alors compléter le graphique avec la densité théorique $f(x)=ae^{−a x}$.

  1. from matplotlib.pyplot import *
  2. from random import *
  3. from math import log,exp
  4. def Expo(a):
  5. return -1/a*log(1-random())
  6. def tirages(n,a):
  7. L=[]
  8. for i in range(n):
  9. L.append(Expo(a))
  10. return L
  11. def graphe(n,a):
  12. L=tirages(n,a)
  13. maxL=int(max(L))
  14. hist(L , bins = 100, range=(0,maxL), normed = True)
  15. x = [float(i)/100 for i in range(100*maxL)]
  16. f = [a*exp(-a*t) for t in x]
  17. plot(x,f)
  18. graphe(10000,0.5)

Télécharger


Enfin l’on peut s’intéresser à la fonction de répartition, que l’on peut simuler avec un histogramme des effectifs cumulés croissants de densité 1 :

  1. def repartition(n,a):
  2. L=tirages(n,a)
  3. maxL=int(max(L))
  4. hist(L , bins = 100, range=(0,maxL), density = True,cumulative=True)
  5. show()

Télécharger

Utilisation :

  1. >>> repartition(10000,0.4)


On peut alors compléter le graphique avec la fonction de répartition théorique $F(x)=1−e^{−ax}$.
Il y a un paramètre dans la commande hist de matplotlib : cumulative=True qui permet d’avoir les Effectifs Cumulés Croissants et avec density=True on a les Fréquences Cumulées Croissantes, donc la fonction de répartition.

  1. def repartition(n,a):
  2. L=tirages(n,a)
  3. maxL=int(max(L))
  4. hist(L , bins = 100, range=(0,maxL), density = True,cumulative=True)
  5. x = [float(i)/100 for i in range(100*maxL)]
  6. f = [1-exp(-a*t) for t in x]
  7. plot(x,f)
  8. show()

Télécharger

Utilisation :

  1. >>> repartition(50000,3)

Demi-vie d’un échantillon de grande taille d’atomes radioactifs.

Décroissance radioactive :
Pour simuler le caractère aléatoire de la désintégration d’un noyau individuel, on peut simuler à l’aide d’un programme Python, si l’on connait la probabilité p de désintégration d’un noyau, on utilise une variable aléatoire qui, si elle est inférieure à p élimine le noyau. Pour travailler sur un échantillon, on fait cela pour n noyaux.

  1. from random import *
  2. def compte_noyaux_restants(n,p):
  3. N=n
  4. for i in range(n):
  5. if random() < p :
  6. N=N-1
  7. return N

Télécharger

Utilisation :

  1. >>> compte_noyaux_restants(10000,0.02)
  2. 9821

Télécharger

On peut alors réitérer cela jusqu’á ce qu’il n’y ait plus de noyaux (En fait on s’arrêtera lorsqu’il en reste moins de 1%) :

  1. def Liste_nb_noyaux_restants(n,p):
  2. L=[]
  3. N=n
  4. while n>0.01*N:
  5. L.append(n)
  6. n=compte_noyaux_restants(n,p)
  7. return L

Télécharger

Utilisation :

  1. >>> Liste_nb_noyaux_restants(10000,0.02)
  2. [10000, 9802, 9606, 9417, 9240, 9069, 8886, ...., 111, 108, 106, 102, 102]

Télécharger

Et la représentation graphique :

  1. import matplotlib. pyplot as plt
  2. def simule_desintegration(n,p):
  3. L=Liste_nb_noyaux_restants(n,p)
  4. plt.plot(range(len(L)),L, linestyle ="none",marker=".")
  5. plt.show()
  6. simule_desintegration(10000,0.03)

Télécharger


Demi-vie :
La demi-vie d’un noyau radioactif, également appelée période radioactive, est la durée nécessaire pour que la moitié des noyaux initialement présents dans un échantillon macroscopique se soit désintégrée. En raison de l’absence de vieillissement, cette demi-vie, caractéristique du noyau considéré, est indépendante de l’instant initial. La demi-vie est, généralement, notée $T_{\frac{1}{2}}$.

  1. def demi_vie(N,p):
  2. n=N
  3. i=0
  4. while n>N/2:
  5. i+=1
  6. n=compte_noyaux_restants(n,p)
  7. return i

Télécharger

Utilisation :

  1. >>> demi_vie(10000,0.03)
  2. 23
  3. >>> demi_vie(1000,0.02)
  4. 36
  5. >>> demi_vie(10000,0.02)
  6. 35

Télécharger

Datation au carbone 14 :
Les organismes vivants contiennent naturellement du carbone 14 ($^{14}C$) provenant de l’interaction des rayons cosmiques avec le carbone présent dans l’atmosphère. La proportion de $^{14}C$ par rapport au $^{12}C$ présent dans un organisme vivant est constante, égale à $1,3\times 10^{-12}$. À la mort d’un organisme, les échanges avec l’atmosphère cessent. Le carbone 14 qu’il contient se désintègre à raison de 12 pour 1000 tous les 100 ans, alors que le $^{12}C$ n’évolue pas.
La demi-vie du $^{14}C$ est d’environ 5800 ans (5 730 ± 40 ans selon Wikipedia) :

  1. >>> demi_vie(100000,12/1000)
  2. 58
  3. >>> demi_vie(100000,12/100000)
  4. 5794

Télécharger

La méthode la plus courante de datation consiste à déterminer la concentration $C_t$ de radiocarbone (c’est-à-dire le rapport $\frac{^{14}C}C_{total}$) d’un échantillon à l’instant $t$ de mesure ; l’âge de l’échantillon est alors donné par la formule : $t − t_0 = \frac{1}{λ}\times ln ⁡\frac{C_0}{C_t}$ où ${C}_{0}$ est la concentration de radiocarbone de l’échantillon à l’instant $t_0$ de la mort de l’organisme d’où provient l’échantillon $\left(C_0 \approx 1.3\times 10^{− 12}\right)$ et λ la constante radioactive du carbone 14$\left(\lambda = \frac {\ln 2}{t_{\frac {1}{2}}} \approx 1,210\cdot 10^{-4} {an} ^{-1}\right)$.
D’où l’algorithme :

  1. from math import log
  2. def datation_C14(p):
  3. C0=1.3*10**(-12)
  4. lamb=1.21*10**(-4)
  5. return round(1/lamb*log(C0/p))

Télécharger

Des archéologues ont trouvé des fragments d’os dans lesquels la proportion de carbone 14 par rapport au carbone 12 est égale à $7\times 10^{-13}$. Estimer l’âge de ces fragments d’os.

  1. >>> datation_C14(7*10**(-13))
  2. 5116

Télécharger

On a découvert dans une grotte en Dordogne un foyer contenant du charbon de bois. À quantité égale, un morceau de bois actuel contient 1,5 fois plus de $^{14}C$ que le charbon de bois prélevé dans la grotte. Estimer une datation de l’occupation de la grotte.

  1. >>> datation_C14(1.3*10**(-12)/1.5)
  2. 3351

Télécharger

Lorsque la teneur en carbone 14 d’un organisme devient inférieure à 0,3 % de sa valeur initiale, on ne peut plus dater précisément à l’aide du carbone 14. Déterminer l’âge à partir duquel un organisme ne peut plus être daté au carbone 14.

  1. >>> datation_C14(0.3/100*1.3*10**(-12))
  2. 48009

Télécharger

Sur des exemples, résolution approchée d’une équation différentielle par la méthode d’Euler.

Résolution par la méthode d’Euler de $y’ = y$, et de $y’ = ay + b$.

En mathématiques, la méthode d’Euler est une procédure numérique pour résoudre par approximation des équations différentielles du premier ordre avec une condition initiale. C’est la plus simple des méthodes de résolution numérique des équations différentielles.
Pour h proche de 0, on a $y(a+h) \approx y(a) + h y’(a)$ .Nous allons utiliser cette approximation affine pour construire pas à pas une fonction vérifiant une équation différentielle du premier ordre et passant par un point donné $(x_0 ; y_0)$.
Soit l’équation différentielle définie par $y’=y$ et les conditions initiales $(x_0 ; y_0)$. En $(x_0 ; y_0)$, on connaît la pente de la tangente à partir de l’équation différentielle, $y_0$. On assimile alors sur l’intervalle $[x_0 ; x_0 + h]$ la fonction à sa tangente. On détermine alors le point $(x_1 ; y_1)$ avec $x_1 = x_0 + h$ et $y_1 = y_0 + hy_0 = y_0(1 + h)$. On recommence le même raisonnement avec le point $(x_1 ; y_1)$. On poursuit en construisant la suite de points $(x_n ; y_n)$ et en assimilant la courbe à une application affine par morceaux.

  1. def Euler(x0,xf,y0,n):
  2. x=x0
  3. y=y0
  4. h=(xf-x0)/float(n)
  5. abscisse=[x0]
  6. ordonnee=[y0]
  7. for i in range(n):
  8. x=x+h
  9. y=y*(1+h)
  10. abscisse.append(x)
  11. ordonnee.append(y)
  12. return ordonnee

Télécharger

Utilisation :

  1. >>> Euler(0,3,1,30)
  2. [1, 1.1, 1.2100000000000002, 1.3310000000000004, 1.4641000000000006, 1.6105100000000008, 1.771561000000001,
  3. 1.9487171000000014, 2.1435888100000016, 2.357947691000002, 2.5937424601000023, 2.853116706110003,
  4. 3.1384283767210035, 3.4522712143931042, 3.797498335832415, 4.177248169415656, 4.594972986357222,
  5. 5.054470284992944, 5.559917313492239, 6.115909044841463, 6.72749994932561, 7.400249944258172, 8.140274938683989,
  6. 8.954302432552389, 9.849732675807628, 10.834705943388391, 11.91817653772723, 13.109994191499954,
  7. 14.420993610649951, 15.863092971714948, 17.449402268886445]

Télécharger

Le graphique ci-dessous est réalisé avec les données ci-dessus et les valeurs exactes.

Soit l’équation différentielle définie par $y’= ay + b$ et les conditions initiales $(x_0 ; y_0)$. En $(x_0 ; y_0)$, on connaît la pente de la tangente à partir de l’équation différentielle, $ay_0 + b$. On assimile alors sur l’intervalle $[x_0 ; x_0 + h]$ la fonction à sa tangente. On détermine alors le point $(x_1 ; y_1)$ avec $x_1 = x_0 + h$ et $y_1 = y_0 + h (ay_0 + b)$. On recommence le même raisonnement avec le point $(x_1 ; y_1)$. On poursuit en construisant la suite de points $(x_n ; y_n)$ et en assimilant la courbe à une application affine par morceaux.

  1. def Euler(x0,xf,y0,n,a,b):
  2. x=x0
  3. y=y0
  4. h=(xf-x0)/float(n)
  5. abscisse=[x0]
  6. ordonnee=[y0]
  7. for i in range(n):
  8. x=x+h
  9. y=y+h*(a*y+b)
  10. abscisse.append(x)
  11. ordonnee.append(y)
  12. return ordonnee

Télécharger

Utilisation :

  1. >>> Euler(0,3,2,30,-0.5,0)
  2. [2, 1.9, 1.805, 1.71475, 1.6290125, 1.547561875, 1.47018378125, 1.3966745921875001, 1.3268408625781252,
  3. 1.2604988194492188, 1.1974738784767578, 1.13760018455292, 1.080720175325274, 1.0266841665590103,
  4. 0.9753499582310597, 0.9265824603195068, 0.8802533373035314, 0.8362406704383548, 0.7944286369164371,
  5. 0.7547072050706152, 0.7169718448170844, 0.6811232525762302, 0.6470670899474187, 0.6147137354500477,
  6. 0.5839780486775453, 0.5547791462436681, 0.5270401889314846, 0.5006881794849104, 0.4756537705106649,
  7. 0.45187108198513165, 0.4292775278858751]

Télécharger


Et l’on peut généraliser le processus :
Soit l’équation différentielle définie par $y’=f(x,y)$ et les conditions initiales $(x_0 ; y_0)$. En $(x_0 ; y_0)$, on connaît la pente de la tangente à partir de l’équation différentielle, $f(x_0 ; y_0)$. On assimile alors sur l’intervalle $[x_0 ; x_0 + h]$ la fonction à sa tangente. On détermine alors le point $(x_1 ; y_1)$ avec $x_1 = x_0 + h$ et $y_1 = y_0 + hf(x_0,y_0)$. On recommence le même raisonnement avec le point $(x_1 ; y_1)$. On poursuit en construisant la suite de points $(x_n ; y_n)$ et en assimilant la courbe à une application affine par morceaux.
La fonction "Euler" ci-dessous prend comme dernier argument une fonction de deux variables :

  1. def Euler(x0,xf,y0,n,f):
  2. x=x0
  3. y=y0
  4. h=(xf-x0)/float(n)
  5. abscisse=[x0]
  6. ordonnee=[y0]
  7. for i in range(n):
  8. x=x+h
  9. y=y+h*f(x,y)
  10. abscisse.append(x)
  11. ordonnee.append(y)
  12. return ordonnee

Télécharger

On peut tester la fonction précédente avec f définie par f(x,y)=-0.5y :

  1. >>> def f(x,y):
  2. ... return -0.5*y
  3. ...
  4. >>> Euler(0,3,2,30,f)
  5. [2, 1.9, 1.805, 1.71475, 1.6290125, 1.547561875, 1.47018378125, 1.3966745921875001, 1.3268408625781252,
  6. 1.2604988194492188, 1.1974738784767578, 1.13760018455292, 1.080720175325274, 1.0266841665590103,
  7. 0.9753499582310597, 0.9265824603195068, 0.8802533373035314, 0.8362406704383548, 0.7944286369164371,
  8. 0.7547072050706152, 0.7169718448170844, 0.6811232525762302, 0.6470670899474187, 0.6147137354500477,
  9. 0.5839780486775453, 0.5547791462436681, 0.5270401889314846, 0.5006881794849104, 0.4756537705106649,
  10. 0.45187108198513165, 0.4292775278858751]

Télécharger

Ce qui donne bien le résultat précédent comme on s’y attendait.

Simulation d’une variable de Bernoulli ou d’un lancer de dé (ou d’une variable uniforme sur un ensemble fini) à partir d’une variable aléatoire de loi uniforme sur [0,1].

Simulation d’une variable de Bernoulli
Une épreuve de Bernoulli a deux issues de probabilités $p$ et $1-p$.
La fonction epreuve() ci-dessous simule une épreuve de Bernoulli :

  1. from random import random
  2. def epreuve(p):
  3. if random()<=p:
  4. return True
  5. else:
  6. return False

Télécharger

La fonction schema() ci-dessous simule un schéma de Bernoulli, elle renvoie le nombre de succés lors de n épreuves de Bernoulli :

  1. def schema(n,p):
  2. nb_succes=0
  3. for i in range(n):
  4. if epreuve(p)==True:
  5. nb_succes+=1
  6. return nb_succes

Télécharger

Pour la simulation de N schémas de Bernoulli avec stockage de la fréquence du succès, on utilisera la fonction experience() ci-dessous :

  1. def experience(N,n,p):
  2. frequences=[]
  3. for i in range(N):
  4. frequences.append(0)
  5. for i in range(N):
  6. nb_succes=schema(n,p)
  7. frequences[i]=nb_succes/n
  8. return frequences

Télécharger

Utilisation :

  1. >>> experience(10,1000,0.4)
  2. [0.4, 0.389, 0.394, 0.406, 0.399, 0.4, 0.401, 0.404, 0.386, 0.415]

Télécharger

Simulation d’un lancer de dé

  1. from random import random
  2. def lancer_de(n):
  3. de=[0]*6
  4. for i in range(n):
  5. valeur=random()
  6. if valeur<=1/6:
  7. de[0]+=1
  8. elif valeur<=2/6:
  9. de[1]+=1
  10. elif valeur<=3/6:
  11. de[2]+=1
  12. elif valeur<=4/6:
  13. de[3]+=1
  14. elif valeur<=5/6:
  15. de[4]+=1
  16. else:
  17. de[5]+=1
  18. return de

Télécharger

Utilisation :

  1. >>> lancer_de(1000)
  2. [177, 177, 143, 170, 178, 155]

Télécharger

Simulation du comportement de la somme de n variables aléatoires indépendantes et de même loi.

Soit n un nombre entier naturel supérieur ou egal à 2. $X_1, X_2, ... , X_n$ sont des variables aléatoires qui suivent la loi uniforme sur [0 ; 1].
Soit Y la variable aléatoire qui prend toutes les valeurs $x_1 + x_2 + ... + x-n$ correspondant à la somme des valeurs prises par $X_1, X_2, ... , X_n$. Y est donc à valeur dans [0 ; n], quelle est la loi suivie par Y ?
On peut utiliser une fonction Python pour simuler cette variable aléatoire :

  1. from random import *
  2. def somme(N,n):
  3. L=[0]*n
  4. for i in range(N):
  5. x=0
  6. for k in range(n):
  7. x=x+random()
  8. x=int(x)
  9. L[x]+=1
  10. for r in range(n):
  11. L[r]=L[r]/1000
  12. return L

Télécharger

Après avoir sommé les n valeurs, on les arrondit afin d’obtenir des effectifs pour n valeurs discrètes. On peut ensuite tracer l’allure de la fonction de répartition de Y :

  1. from random import *
  2. import matplotlib.pyplot as pyplot
  3. def somme(N,n):
  4. L=[0]*n
  5. for i in range(N):
  6. x=0
  7. for k in range(n):
  8. x=x+random()
  9. x=round(x)
  10. L[x]+=1
  11. for r in range(n):
  12. L[r]=L[r]/1000
  13. A=[k for k in range (n)]
  14. pyplot.bar(A,L,width=0.5,color="blue")
  15. pyplot.show()
  16. return L

Télécharger

Utilisation :

  1. >>> somme(1000,16)
  2. [0, 0, 0, 0, 3, 43, 142, 349, 279, 147, 35, 2, 0, 0, 0, 0]
  3. [0.0, 0.0, 0.0, 0.0, 0.003, 0.043, 0.142, 0.349, 0.279, 0.147, 0.035, 0.002, 0.0, 0.0, 0.0, 0.0]

Télécharger


Version 2 en utilisant la fonction hist de matplotlib :

  1. from random import *
  2. import matplotlib.pyplot as pyplot
  3. def somme2(N,n):
  4. L=[]
  5. for i in range(N):
  6. x=0
  7. for k in range(n):
  8. x=x+random()
  9. L.append(x)
  10. pyplot.hist(L,bins=100,range=(1,n+1))
  11. pyplot.show()

Télécharger

Utyilisation :

  1. >>> somme2(1000,16)


notes

[2La fonction hist (si on a écrit from matplotlib.pyplot import * en pré-ambule) trace un histogramme d’une liste donnée en argument (compte le nombre d’éléments de la liste dans différents intervalles). On peut rajouter quelques arguments supplémentaires : par exemple hist(liste , bins =10, range=(0,4), normed = True) où bins donne le nombre de ’bâtons’ de l’histogramme, range impose l’intervalle sur lequel on fait l’histogramme, normed=True renormalise l’histogramme pour que l’aire totale des bâtons soit égale à 1.

Réagir à cet article
Vous souhaitez compléter cet article pour un numéro futur, réagir à son contenu, demander des précisions à l'auteur ou au comité de rédaction...
À lire aussi ici
MathémaTICE est un projet
en collaboration avec
Suivre la vie du site Flux RSS 2.0  |  Espace de rédaction  |  Nous contacter  |  Site réalisé avec: SPIP  |  N° ISSN 2109-9197