Mathématice, intégration des Tice dans l'enseignement des mathématiques  
Sommaire > N°71 - Septembre 2020 - En cours d’élaboration > Les algorithmes du programme de l’option de (...)

Les algorithmes du programme de l’option de Terminale "Mathématiques expertes" (2020).
Moteur de recherche

On peut lire dans le programme officiel de cette option :

L’utilisation de logiciels (calculatrice ou ordinateur), d’outils de visualisation et de représentation, de calcul (numérique ou formel), de simulation, de programmation développe la possibilité d’expérimenter, favorise l’interaction entre l’observation et la démonstration,et change profondément la nature de l’enseignement. L’utilisation régulière de ces outils peut intervenir selon trois modalités :

  • par le professeur, en classe, avec un dispositif de visualisation collective adapté ;
  • par les élèves, sous forme de travaux pratiques de mathématiques en classe, à l’occasion de la résolution d’exercices ou de problèmes ;
  • dans le cadre du travail personnel des élèves hors du temps de classe (par exemple au CDI ou à un autre point d’accès au réseau local).

Des "Exemples d’algorithmes" sont proposés :

Algorithme d’Euclide de calcul du PGCD de deux nombres et calcul d’un couple de Bézout.

Algorithme d’Euclide de calcul du PGCD de deux nombres :
On peut commencer par l’algorithme des différences :

  1. def pgcd(a,b):
  2. if b==0:
  3. return a
  4. else:
  5. while a!=b:
  6. if a>=b:
  7. a -= b
  8. else:
  9. b -= a
  10. return a

Télécharger

Puis l’algorithme des divisions successives :

  1. def pgcd(a,b):
  2. if b==0:
  3. return a
  4. else:
  5. while b!=0:
  6. a,b=b,a%b
  7. return a

Télécharger

Le même en version récursive :

  1. def pgcd(a,b):
  2. if b==0:
  3. return a
  4. else:
  5. return pgcd(b,a%b)

Télécharger

Calcul d’un couple de Bézout :

  1. def bezout(a,b):
  2. (u0,v0,u1,v1)=(1,0,0,1)
  3. while b:
  4. (q,r)=divmod(a,b)
  5. (a,b)=(b,r)
  6. (u0,v0,u1,v1)=(u1,v1,u0-q*u1,v0-q*v1)
  7. return (u0,v0)

Télécharger

Utilisation :

  1. >>> bezout(412,326)
  2. (-72, 91)

Télécharger

Version récursive :

  1. def bezout(a, b):
  2. if b == 0:
  3. return (1,0)
  4. else:
  5. (u,v) = bezout(b,a%b)
  6. return (v, u-(a//b)*v)

Télécharger

Crible d’Ératosthène.

Le crible d’Ératosthène est une méthode qui fournit la liste des nombres premiers inférieurs à une valeur fixée n :
- on forme la liste des entiers de 2 à n ;
- on retient comme « nombre premier » le premier nombre de la liste non encore barré (le premier dans ce cas est 2) ;
- on barre tous les entiers multiples du nombre retenu à l’étape précédente, en commençant par son carré (puisque 2 × i, 3 × i, … , (i – 1) × i ont déjà été barrés en tant que multiples de 2, 3, ...) ;
- on répète les deux dernières opérations (c’est-à-dire : on retient le prochain nombre non barré et on barre ses multiples) ;
- dès qu’on en est à chercher les multiples des nombres excédant la racine carrée de n, on termine l’algorithme.
Ainsi les nombres premiers inférieurs à n sont les nombres qui restent non barrés à la fin du processus.
Pour programmer cela en Python, on crée une liste [0,0,2,3,4,...,n] et au lieu de barrer les multiples, on va les remplacer par des 0, à la fin on retourne la liste sans les 0 :

  1. def eratosthene(n):
  2. """retourne le tableau des nombres premiers <= n (crible d'eratosthene)"""
  3. n += 1
  4. tableau = [0,0]+[i for i in range(2,n)]
  5. for i in range(2,n):
  6. if tableau[i]!= 0:
  7. #c'est un nombre 1er : on le garde, mais on supprime ses multiples
  8. for j in range(i*2,n,i):
  9. tableau[j] = 0
  10. return [p for p in tableau if p!=0]

Télécharger

On peut améliorer ce code :

  1. def eratosthene2(n):
  2. """retourne la liste des nombres premiers <= n (crible d'eratosthene)"""
  3. if n<2:
  4. return []
  5. n += 1
  6. tableau = [False,False] + [True]*(n-2)
  7. tableau[2::2] = [False]*((n-2)//2 + n%2) # sup. des nb pairs
  8. premiers = [2] # initialisation du tableau des nb 1ers (2 est 1er)
  9. racine = int(n**0.5)
  10. racine = racine + [1,0][racine%2] # pour que racine soit impair
  11. for i in range(3, racine+1, 2):
  12. if tableau[i]:
  13. # on enregistre le nouveau nb 1er
  14. premiers.append(i)
  15. # on élimine i et ses multiples
  16. tableau[i::i] = [False]*((n-i)//i + int((n-i)%i>0))
  17. for i in range(racine, n, 2):
  18. if tableau[i]:
  19. # on enregistre le nouveau nb 1er (pas de multiples à supprimer)
  20. premiers.append(i)
  21. return premiers

Télécharger

Décomposition en facteurs (...)

Décomposition en facteurs premiers

  1. def facteurs(n):
  2. F=[]
  3. x=n
  4. i=2
  5. while i<n//2+1:
  6. a,b=divmod(x,i)
  7. if b==0:
  8. F.append(i)
  9. x=a
  10. else:
  11. i+=1
  12. if F==[]:
  13. F.append(n)
  14. return F

Télécharger

On peut accélérer les calculs en arrangeant un peu le code :

  1. from math import sqrt
  2. def facteurs(n):
  3. """facteurs(n): décomposition d'un nombre entier n en facteurs premiers"""
  4. F = []
  5. if n==1:
  6. return F
  7. # recherche de tous les facteurs 2 s'il y en a
  8. while n>=2:
  9. x,r = divmod(n,2)
  10. if r!=0:
  11. break
  12. F.append(2)
  13. n = x
  14. # recherche des facteurs 1er >2
  15. i=3
  16. rn = int(sqrt(n))
  17. while i<=n:
  18. if i>rn:
  19. F.append(n)
  20. break
  21. x,r = divmod(n,i)
  22. if r==0:
  23. F.append(i)
  24. n=x
  25. rn = int(sqrt(n))
  26. else:
  27. i += 2
  28. return F

Télécharger

Algorithme PageRank

Le PageRank ou PR est l’algorithme d’analyse des liens concourant au système de classement des pages Web utilisé par le moteur de recherche Google. Il mesure quantitativement la popularité d’une page web. Le PageRank n’est qu’un indicateur parmi d’autres dans l’algorithme qui permet de classer les pages du Web dans les résultats de recherche de Google. Ce système a été inventé par Larry Page, cofondateur de Google. Ce mot est une marque déposée.
Le principe de base est d’attribuer à chaque page une valeur (ou score) proportionnelle au nombre de fois que passerait par cette page un utilisateur parcourant le graphe du Web en cliquant aléatoirement, sur un des liens apparaissant sur chaque page. Ainsi, une page a un PageRank d’autant plus important qu’est grande la somme des PageRanks des pages qui pointent vers elle (elle comprise, s’il y a des liens internes). Le PageRank est une mesure de centralité sur le réseau du web.
Plus formellement, le déplacement de l’utilisateur est une marche aléatoire sur le graphe du Web, c’est-à-dire le graphe orienté dont les sommets représentent les pages du Web et les arcs les hyperliens. En supposant que l’utilisateur choisisse chaque lien indépendamment des pages précédemment visitées (le réalisme d’une telle hypothèse pouvant être discuté).
Supposons un site qui possède une page d’accueil et deux pages (page 1 et page 2), on peut naviguer comme l’on veut dans le site sauf aller de la page 2 à la page 1 :

Graphe réalisé en ligne sur le site Graphviz. Algorithme inspiré de SNT Numworks.
D’où l’algorithme suivant :

  1. from random import *
  2. def PR(n):
  3. nAccueil,n1,n2=0,0,0
  4. start=randint(0,2)
  5. if start==0:
  6. page="Accueil"
  7. nAccueil=1
  8. elif start==1:
  9. page="page 1"
  10. n1=1
  11. else:
  12. page="page 2"
  13. n2=1
  14. for i in range(1,n):
  15. if page=="Accueil":
  16. alea=randint(0,1)
  17. if alea==0:
  18. page="page 1"
  19. n1=n1+1
  20. else:
  21. page="page 2"
  22. n2=n2+1
  23. elif page=="page 1":
  24. page="Accueil"
  25. nAccueil=nAccueil+1
  26. else:
  27. alea=randint(0,1)
  28. if alea==0:
  29. page="Accueil"
  30. nAccueil=nAccueil+1
  31. else:
  32. page="page 1"
  33. n1=n1+1
  34. fAccueil,f1,f2=round(nAccueil/n,2),round(n1/n,2),round(n2/n,2)
  35. return fAccueil,f1,f2

Télécharger

Utilisation :

  1. >>> PR(100)
  2. (0.46, 0.33, 0.21)
  3. >>> PR(1000)
  4. (0.45, 0.34, 0.21)
  5. >>> PR(100000)
  6. (0.44, 0.33, 0.22)

Télécharger

En dehors de ces exemples répertoriés, il est bien sûr possible de proposer des programmes Python pour les différents chapitres.

Les nombres complexes

Python permet le calcul avec des nombres complexes (Attention cependant Python utilise la lettre j là où nous utilisons la lettre i ...) :

  1. >>> z1=complex(-2,3)
  2. >>> z1
  3. (-2+3j)
  4. >>> z1.real
  5. -2.0
  6. >>> z1.imag
  7. 3.0
  8. >>> z1.conjugate()
  9. (-2-3j)
  10. >>> z2=complex(1,-1)
  11. >>> z1+z2
  12. (-1+2j)
  13. >>> z1-z2
  14. (-3+4j)
  15. >>> z1*z2
  16. (1+5j)
  17. >>> z1/z2
  18. (-2.5+0.5j)
  19. >>> z2**4
  20. (-4-0j)
  21. >>> 1/complex(1,1)
  22. (0.5-0.5j)
  23. >>> abs(complex(-1,0))
  24. 1.0
  25. >>> abs(complex(1,1))
  26. 1.4142135623730951

Télécharger

abs(z) renvoie le module du nombre complexe z.
Et en utilisant le module cmath (Fonctions mathématiques pour nombres complexes) :

  1. from cmath import *
  2. >>> phase(z1)
  3. 2.158798930342464
  4. >>> phase(complex(-1,0))
  5. 3.141592653589793
  6. >>> phase(complex(0,1))
  7. 1.5707963267948966
  8. >>> polar(complex(10,0))
  9. (10.0, 0.0)
  10. >>> polar(complex(0,-1))
  11. (1.0, -1.5707963267948966)
  12. >>> polar(complex(-1,0))
  13. (1.0, 3.141592653589793)
  14. >>> rect(2,pi/2)
  15. (1.2246467991473532e-16+2j)
  16. >>> rect(2,pi/3)
  17. (1.0000000000000002+1.7320508075688772j)

Télécharger

phase(z) renvoie l’argument du nombre complexe z appartenant à l’intervalle $]-\pi ; \pi]$.
polar(z) renvoie une paire (r, phi) où r est le module de z et phi est l’argument de z appartenant à l’intervalle $]-\pi ; \pi]$. polar(z) est équivalent à (abs(z), phase(z)).
rect(r,phi) renvoie le nombre complexe z dont les coordonnées polaires sont r et phi. Équivalent à r*(cos(phi) +sin(phi)*1j).
Argument de mesure principale
On peut reprogrammer cette fonction phase en utilisant la fonction $arccos(x)$, définie sur [-1 ; 1] à valeurs dans $[0 ; \pi]$ :

  1. def arg(z):
  2. r=abs(z)
  3. c=z.real/r
  4. if z.imag>=0:
  5. a=math.acos(c)
  6. else:
  7. a=-math.acos(c)
  8. return a

Télécharger

Utilisation et vérification :

  1. >>> arg(complex(1,-1))
  2. -0.7853981633974484
  3. >>> arg(complex(1,1))
  4. 0.7853981633974484
  5. >>> arg(complex(-1,1))
  6. 2.356194490192345
  7. >>> arg(complex(-1,-1))
  8. -2.356194490192345
  9. >>> phase(complex(1,-1))
  10. -0.7853981633974483
  11. >>> phase(complex(1,1))
  12. 0.7853981633974483
  13. >>> phase(complex(-1,-1))
  14. -2.356194490192345
  15. >>> phase(complex(-1,1))
  16. 2.356194490192345

Télécharger

Racine complexe
Pour utiliser ces fonctionnalités, on peut demander aux élèves de proposer une fonction Python qui prend en argument un nombre complexe et qui renvoie ses deux racines carrées sous forme de nombres complexes :

  1. def sqrt_complex(c):
  2. r=sqrt(abs(c))
  3. phi=phase(c)/2
  4. x=r*cos(phi)
  5. y=r*sin(phi)
  6. return complex(-x,-y),complex(x,y)

Télécharger

Utilisation :

  1. >>> sqrt_complex(complex(0,1))
  2. ((-0.7071067811865476-0.7071067811865475j), (0.7071067811865476+0.7071067811865475j))
  3. >>> sqrt_complex(complex(1,1))
  4. ((-1.0986841134678098-0.45508986056222733j), (1.0986841134678098+0.45508986056222733j))

Télécharger

En effet, $\sqrt{a + ib} = \sqrt{re^{i\theta}} = \sqrt{r}e^{i\frac \theta 2}$, d’où le script ci-dessus.
Il est également possible d’avoir une solution et un script purement algébrique en remarquant que si $(x + iy)^2 = a+ ib$ alors $x^2 - y^2 = a$ et $2xy = b$, de plus $x^2 + y^2 = \sqrt{a^2 + b^2}$, alors $x^2 = \frac{\sqrt{a^2 + b^2} + a} 2$ et $y^2 = \frac{\sqrt{a^2 + b^2} - a} 2$.
Et puisque $2xy = b$, si b>0 : x et y sont de même signe, sinon ils sont de signes contraires, d’où le script suivant :

  1. def sqrt_complex2(c):
  2. a=c.real
  3. b=c.imag
  4. module=abs(c)
  5. x=sqrt((a+module)/2)
  6. y=sqrt((module-a)/2)
  7. if b<0:
  8. y=-y
  9. return complex(-x,-y),complex(x,y)

Télécharger

Solutions dans C d’une équation du second degré à coefficients réels

  1. def polynome_degre2(a,b,c):
  2. delta=b**2-4*a*c
  3. if delta>0:
  4. print("1")
  5. return ["Deux solutions réelles :",(-b-math.sqrt(delta))/(2*a),(-b+math.sqrt(delta))/(2*a)]
  6. elif b==0:
  7. print("2")
  8. return ["Une solution réelle :",-b/(2*a)]
  9. else:
  10. print("3")
  11. return ["Deux solutions complexes :",complex(-b/(2*a),-math.sqrt(-delta)/(2*a)),complex(-b/(2*a),math.sqrt(-delta)/(2*a))]

Télécharger

Utilisation :

  1. >>> polynome_degre2(1,1,0)
  2. ['Deux solutions réelles :', -1.0, 0.0]
  3. >>> polynome_degre2(25,-4,1/4)
  4. ['Deux solutions complexes :', (0.08-0.06j), (0.08+0.06j)]
  5. >>> polynome_degre2(1,-14,49)
  6. ['Une solution réelle :', 7.0]

Télécharger

Algorithme de Horner
Lorsqu’on cherche à évaluer un polynôme $P(x)= a_n x^n + a_{n−1}x^{n−1} + ⋯ + a_0$, il n’est clairement pas optimal de calculer chaque $a_k x^k$, puis de les ajouter. En effet, pour calculer $a_n x^n$, il est intéressant de se souvenir que l’on a déjà calculé $a_{n−1}x^{n−1}$, et donc $x^{n−1}$.
Le mathématicien anglais Horner a mis au point une méthode efficace utilisant cette remarque. Voici comment elle fonctionne sur un exemple simple : pour calculer $5x^4−4x^3+3x^2−2x+1$, on peut remarquer que cette quantité est égale à $x(x(x(5x−4)+3)−2)+1$. On calcule alors successivement $a=5, b=xa−4=5x−4, c = xb + 3= x(5x − 4) + 3, d = xc − 2$, et enfin $e = xd + 1$ qui est le résultat souhaité. Plus généralement, pour un polynôme $P(x)=a_n x^n + a_{n−1}x^{n−1}+⋯+a_1 x+a_0$, on écrit $P(x)=x(x(⋯x(a_n x + a_{n−1}) + a_{n−2})⋯) + a_0$.
Le calcul nécessitera simplement $n$ multiplications et $n$ additions contre $\frac{n(n - 1)}2$ multiplications et $n$ additions.
Voici une fonction écrite sous Python qui prend en entrée la liste des coefficients d’un polynôme $P$ et un complexe $z$, et qui retourne $P(z)$ grâce à la méthode de Horner.

  1. def Horner(L,a):
  2. n=len(L)
  3. b=L[0]
  4. for i in range(1,n):
  5. b=b*a+L[i]
  6. return b

Télécharger

Utilisation (avec des complexes ou des réels) :

  1. >>> Horner([2,-5,1,3],10)
  2. 1513
  3. >>> Horner([2,-5,1,3],1j)
  4. (8-1j)
  5. >>> Horner([2,-5,1,3],complex(1,1))
  6. -5j

Télécharger

Utilisation des racines n-ièmes de l’unité pour tracer un polygone régulier
Si n désigne un entier naturel supérieur ou égal à 3, on peut écrire une fonction Python de paramètre n :

  1. from pylab import *
  2. from cmath import *
  3. def polygone(n):
  4. for k in range(n):
  5. z=2*exp(1j*2*k*pi/n)
  6. plot(z.real,z.imag,"r.",markersize=12)
  7. show()

Télécharger

Utilisation :

  1. >>> polygone(10)


Alignement de points
A, B et C sont alignés si et seulement si $\frac{c-a}{b-a}$ est un nombre réel, soit A, B et C sont alignés si et seulement si $arg\left(\frac{c-a}{b-a}\right) \equiv 0 [{\pi}]$, d’où la fonction Python :

  1. from cmath import *
  2. def Points_alignes(a,b,c):
  3. z=(c-a)/(b-a)
  4. if phase(z)==0 or phase(z)==pi:
  5. return True
  6. else:
  7. return False

Télécharger

Utilisation :

  1. >>> Points_alignes(3+5j,-4+4j,-11+3j)
  2. True
  3. >>> Points_alignes(1j,8,-7+2j)
  4. False
  5. >>> Points_alignes(1-1j,5+1j,-1-2j)
  6. True

Télécharger

Droites orthogonales
(AB) et (AC) sont orthogonales si et seulement si $arg\left(\frac{c-a}{b-a}\right) \equiv \frac \pi 2 [{\pi}]$, d’où la fonction Python :

  1. from cmath import *
  2. def Droites_orthogonales(a,b,c):
  3. z=(c-a)/(b-a)
  4. if phase(z)==pi/2 or phase(z)==-pi/2:
  5. return True
  6. else:
  7. return False

Télécharger

>>> Droites_orthogonales(1-1j,11-11j,4+2j)
True
>>> Droites_orthogonales(8-3j,10+6j,6+5j)
False
>>> Droites_orthogonales(6+5j,8-3j,10+6j)
True
Ensemble de Mandelbrot
En mathématiques, l’ensemble de Mandelbrot est une fractale définie comme l’ensemble des points c du plan complexe pour lesquels la suite de nombres complexes définie par récurrence par ${\displaystyle {\begin{cases}z_{0}=0\\z_{n+1}=z_{n}^{2}+c\end{cases}}} $ est bornée.

  1. import matplotlib.pyplot as plt
  2. from random import random
  3. def mandelbrot(n):
  4. for i in range(n):
  5. c=complex(4*random()-2,4*random()-2)
  6. z=0
  7. k=0
  8. while abs(z)<2 and k<50:
  9. z=z**2+c
  10. k=k+1
  11. if k==50:
  12. plt.plot(c.real,c.imag,"rx",ms=0.5)
  13. plt.xlim=(-2,2)
  14. plt.ylim=(-2,2)
  15. ax=plt.gca()
  16. ax.spines["bottom"].set_position(("data",0))
  17. ax.spines["left"].set_position(("data",0))
  18. plt.show()

Télécharger

Utilisation :

  1. >>> mandelbrot(100000)

  1. >>> mandelbrot(500000)

  1. >>> mandelbrot(1000000)


Ensemble de Julia
Un ensemble de Julia $J_c$ est une fractale qui est paramétrée par un nombre complexe $c = a + ib$. Pour déterminer si un point $M(x ; y)$ du plan appartient à l’ensemble de Julia, on définit une suite de couples de nombres $(x_n ; y_n)$ par $(x_0 ; y_0) = (x ; y)$ et par la relation de récurrence : ${\displaystyle {\begin{cases}x_{n+1}=x_n^2 − y_n^2+a\\y_{n+1} = 2x_n \times y_n+b\end{cases}}} $.
Par exemple si $(x_0 ; y_0) = (3 ; 4)$ et si $(a ; b) = (2 ; 1)$ alors :
$x_1 = 3^2 − 4^2 + 2 = −5$ et $y_1 = 2x_0\times y_0 + 1 = 2\times 3\times 2 + 1 = 13$
$x_2 = (−5)^2 − 13^2 + 2 = −142$ et $y_2 = 2\times 1\times y_1 + 1 = 2\times(−5)\times 13 + 1 = −129$
Le point $M(x ; y)$ n’appartient pas à l’ensemble de Julia $J_c$ si le nombre $x_n^2 + y_n^2$ (le carré de la distance du point de coordonnées $(x_n ; y_n)$ à l’origine) tend vers l’infini lorsque $n$ tend vers l’infini. Une propriété permet d’affirmer que si $x_n^2 + y_n^2$ dépasse $2^2$ à partir d’un certain $n$ alors la suite $(x_n^2 + y_n^2)$ diverge vers l’infini et le point initial $M(x ; y)$ n’appartient pas à $J_c$.

  1. from PIL import Image
  2. def Julia(c,taille):
  3. a=c.real
  4. b=c.imag
  5. xmin = -1.25
  6. xmax = 1.25
  7. ymin = -1.25
  8. ymax = 1.25
  9. iterationmax = 200
  10. im = Image.new("RGB",(taille,taille),(255,255,255))
  11. pixels = im.load()
  12. for line in range(taille):
  13. for col in range(taille):
  14. i = 1
  15. x = xmin+col*(xmax-xmin)/taille
  16. y = ymax-line*(ymax-ymin)/taille
  17. while i<=iterationmax and (x**2+y**2)<=4:
  18. stock = x
  19. x = x**2-y**2+a
  20. y = 2*stock*y+b
  21. i += 1
  22. if i>iterationmax and (x**2+y**2)<=4:
  23. pixels[col,line] = (0,0,0)
  24. else:
  25. pixels[col,line] = ((4*i)%256,2*i,(6*i)%256)
  26. im.save("julia1.png")
  27. im.show()

Télécharger

Utilisation (Pensez à changer le nom du fichier sauvegardé si vous ne voulez pas perdre l’image précédente) :

  1. >>> Julia(complex(0.284, 0.0122),400)

  1. >>> Julia(complex(-0.1011,0.9563),400)

  1. Julia(complex(-0.1528,1.0397),400)

  1. Julia(complex(-0.06783611264225832,0.6617460391250546),400)

  1. Julia(complex(-0.577,0.478),800)


Le lapin de Douady
Le lapin de Douady est l’ensemble de Julia obtenu à l’aide du nombre complexe solution de l’équation c³ + 2c² + c + 1 = 0, c’est-à-dire -0.12256 ± 0.74486i. Il a été nommé en hommage à Adrien Douady, mathématicien français mort en 2006.

  1. >>> Julia(complex(-0.12256,0.74486),400)

Arithmétique

Avec le calcul d’un couple de Bézout, on peut alors rajouter la Résolution de l’équation diophantienne ax + by = c :
Une équation diophantienne, en mathématiques, est une équation polynomiale à une ou plusieurs inconnues dont les solutions sont cherchées parmi les nombres entiers relatifs, les coefficients étant eux-mêmes également entiers.

  1. def equation_diophantienne(a,b,c):
  2. p=pgcd(a,b)
  3. x0=c*bezout(a,b)[0]/p
  4. y0=c*bezout(a,b)[1]/p
  5. return [(x0,y0),(b/p,"k +",x0,";"-a/p,"k +",y0)]

Télécharger

Utilisation

  1. >>> equation_diophantienne(3,2,12)
  2. [(12.0, -12.0), (2.0, 'k +', 12.0, ';', -3.0, 'k +', -12.0)]

Télécharger

Traduction : L’équation 3x + 2y = 12 a pour solution particulière (12 ; -12) et pour solution générale : (2k + 12 ; -3k - 12) pour tout k entier relatif.
Inverse modulaire
Dans de nombreuses applications cryptographiques, il est nécessaire de trouver l’inverse modulaire d’un nombre. Soit $0 < a < m$, tel que $a$ et $m$ sont des entiers. L’inverse modulaire de $a$ est l’unique entier $n$ avec $0 < n < m$, tel que le reste dans la division euclidienne de $a\times n$ par $m$ soit égal á 1.
Par exemple, $4\times 13 = 52 = 17 \times 3 + 1$. Alors le reste de la division de 52 par 17 est 1. Ainsi, 13 est l’inverse de 4 modulo 17 (On dira aussi 4 est l’inverse de 13 modulo 17).
Si les nombres $a$ et $m$ ne sont pas premiers entre eux, l’inverse modulaire n’existe pas, d’où le script suivant :

  1. def inverse_modulaire(a,m):
  2. if pgcd(a,m)>1:
  3. return None
  4. else:
  5. return m+bezout(a,m)[0]

Télécharger

Utilisation :

  1. >>> inverse_modulaire(6,31)
  2. 26
  3. >>> inverse_modulaire(6,7)
  4. 6
  5. >>> inverse_modulaire(4,8)
  6. >>> inverse_modulaire(13,17)
  7. 21
  8. >>> inverse_modulaire(39,53)
  9. 34

Télécharger

On peut aussi déterminer cet inverse modulaire sans Bézout :

  1. def inverse_modulaire(x,m):
  2. for n in range(m):
  3. if (x*n)%m == 1:
  4. return n
  5. elif n == m-1:
  6. return "Null"
  7. else:
  8. continue

Télécharger

Diviseurs
Si l’on veut afficher la liste des diviseurs d’un nombre dans $\mathbb{N}$ :

  1. def diviseurs(n):
  2. L=[]
  3. for i in range(1,n+1):
  4. if n%i==0:
  5. L.append(i)
  6. return L

Télécharger

Utilisation :

  1. >>> diviseurs(4567)
  2. [1, 4567]
  3. >>> diviseurs(4568)
  4. [1, 2, 4, 8, 571, 1142, 2284, 4568]
  5. >>> diviseurs(60)
  6. [1, 2, 3, 4, 5, 6, 10, 12, 15, 20, 30, 60]
  7. >>> diviseurs(36)
  8. [1, 2, 3, 4, 6, 9, 12, 18, 36]

Télécharger

Si l’on veut afficher la liste des diviseurs d’un nombre dans $\mathbb{Z}$ :

  1. def diviseursZ(n):
  2. if n<0:
  3. M=diviseurs(-n)
  4. else:
  5. M=diviseurs(n)
  6. L=[]
  7. for x in M:
  8. L.append(-x)
  9. L=L+M
  10. L.sort()
  11. return L

Télécharger

Utilisation :

  1. >>> diviseursZ(72)
  2. [-72, -36, -24, -18, -12, -9, -8, -6, -4, -3, -2, -1, 1, 2, 3, 4, 6, 8, 9, 12, 18, 24, 36, 72]
  3. >>> diviseursZ(-72)
  4. [-72, -36, -24, -18, -12, -9, -8, -6, -4, -3, -2, -1, 1, 2, 3, 4, 6, 8, 9, 12, 18, 24, 36, 72]

Télécharger

Nombre premier
On peut se servir de ce qui précède pour obtenir une fonction qui valide ou non le fait d’être un nombre premier :

  1. def est_premier(n):
  2. L=diviseurs(n)
  3. if len(L)==2:
  4. return True
  5. else:
  6. return False

Télécharger

Utilisation :

  1. >>> est_premier(1)
  2. False
  3. >>> est_premier(2)
  4. True
  5. >>> est_premier(0)
  6. False
  7. >>> est_premier(6)
  8. False
  9. >>> est_premier(19)
  10. True

Télécharger

Liste des n premiers nombres premiers

  1. def liste_premiers(n):
  2. i=0
  3. j=2
  4. premiers=[]
  5. while i<n:
  6. if est_premier(j):
  7. premiers.append(j)
  8. j=j+1
  9. i=i+1
  10. else:
  11. j=j+1
  12. return premiers

Télécharger

Utilisation :

  1. >>> liste_premiers(25)
  2. [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

Télécharger

Densité des nombres premiers
On peut représenter graphiquement tous les nombres entiers de 1 jusqu’a N=m*n et ne faire apparaître que les nombres premiers :

  1. import matplotlib.pyplot as plt
  2. def graphe(m,n):
  3. for colonne in range(1,m+1):
  4. for ligne in range(1,n+1):
  5. nombre =colonne+(ligne-1)*m
  6. if est_premier(nombre)==True:
  7. plt.plot(ligne,colonne,"ko")
  8. plt.xticks([],[])
  9. plt.yticks([],[])
  10. plt.show()

Télécharger

Utilisation :

  1. >>> graphe(5,5)


Les nombres sont ainsi représentés :
5 X X X X
X X X 19 X
3 X 13 X 23
2 7 X 17 X
X X 11 X X
En changeant plt.plot(ligne,colonne,"ko") par plt.plot(ligne,colonne,"k.") :

  1. >>> graphe(100,100)


Nombres amicaux
On appelle diviseur strict de $n$ un nombre entier naturel qui divise $n$, strictement inférieur à $n$ :

  1. def diviseurs_stricts(n):
  2. L=[]
  3. for i in range(1,n//2+1):
  4. if n%i==0:
  5. L.append(i)
  6. return L

Télécharger

Utilisation :

  1. >>> diviseurs_stricts(220)
  2. [1, 2, 4, 5, 10, 11, 20, 22, 44, 55, 110]
  3. >>> diviseurs_stricts(284)
  4. [1, 2, 4, 71, 142]

Télécharger

On dit que les nombres entiers naturels $n$ et $p$, distincts, sont amicaux lorsqu’ils sont respectivement égaux à la somme des diviseurs stricts de l’autre :

  1. def amicaux(n,p):
  2. if sum(diviseurs_stricts(n))==p and sum(diviseurs_stricts(p))==n:
  3. return True
  4. else:
  5. return False

Télécharger

Utilisation :

  1. >>> amicaux(36,48)
  2. False
  3. >>> amicaux(220,284)
  4. True

Télécharger

On peut alors chercher tous les amis inférieurs à n :

  1. def cherche_amis(n):
  2. amis=[]
  3. for i in range(2,n):
  4. for j in range(i+1,n):
  5. if amicaux(i,j)==True:
  6. amis.append((i,j))
  7. return amis

Télécharger

Et ils ne sont pas nombreux :

  1. >>> cherche_amis(1000)
  2. [(220, 284)]
  3. >>> cherche_amis(2000)
  4. [(220, 284), (1184, 1210)]

Télécharger

Nombres parfaits
On dit qu’un nombre est parfait lorsque la somme de ses diviseurs stricts est égale à lui-même, autrement dit lorsqu’il est ami avec lui-même, d’où le script suivant :

  1. def est_parfait(n):
  2. if sum(diviseurs_stricts(n))==n:
  3. return True
  4. else:
  5. return False

Télécharger

Utilisation :

  1. >>> est_parfait(6)
  2. True
  3. >>> est_parfait(28)
  4. True
  5. >>> est_parfait(36)
  6. False

Télécharger

On peut alors chercher tous les parfaits inférieurs à n :

  1. def cherche_parfait(n):
  2. parfaits=[]
  3. for i in range(2,n):
  4. if est_parfait(i)==True:
  5. parfaits.append(i)
  6. return parfaits

Télécharger

Utilisation :

  1. >>> cherche_parfait(500)
  2. [6, 28, 496]
  3. >>> cherche_parfait(5000)
  4. [6, 28, 496]
  5. >>> cherche_parfait(50000)
  6. [6, 28, 496, 8128]

Télécharger

Tableau de congruences
En se servant du script ci-dessous, l’on peut, sans trop de difficulté obtenir un tableau de congruences :

  1. print("{:22}".format("n ≡ ... [4] |"),end="")
  2. for n in range(4):
  3. print("{:2d}".format(n%4),end=" | ")
  4. print()
  5. print("{:22}".format("2n ≡ ... [4] | "),end="")
  6. for n in range(4):
  7. print("{:2d}".format(2*n %4),end=" | ")
  8. print()
  9. print("{:22}".format("n² + 1 ≡ ... [4] | "),end="")
  10. for n in range(4):
  11. print("{:2d}".format((n**2+1)%4),end=" | ")
  12. print()
  13. print("{:22}".format("2n(n² + 1) ≡ ... [4]| "),end="")
  14. for n in range(4):
  15. print("{:2d}".format((2*n*(n**2+1))%4),end=" | ")

Télécharger

Utilisation :

  1. n ≡ ... [4] | 0 | 1 | 2 | 3 |
  2. 2n ≡ ... [4] | 0 | 2 | 0 | 2 |
  3. n² + 1 ≡ ... [4] | 1 | 2 | 1 | 2 |
  4. 2n(n² + 1) ≡ ... [4]| 0 | 0 | 0 | 0 |

Télécharger

Décomposition de la fraction $\frac{n}{d}$ en fraction continue

  1. def frac_continue(n,d):
  2. solution=[]
  3. while d:
  4. a=n//d
  5. solution.append(a)
  6. n,d=d,n-a*d
  7. return solution

Télécharger

Utilisation :

  1. >>> frac_continue(314,100)
  2. [3, 7, 7]
  3. >>> frac_continue(31415,10000)
  4. [3, 7, 14, 1, 8, 2]

Télécharger

Ce que l’on peut traduire par :
$3,14 = \frac{314}{100} = 3 + \frac{1}{7 + \frac{1}{7}}$ et $3,1415 = \frac{31415}{10000} = 3 + \frac{1}{7 + \frac{1}{14 + \frac{1}{1+\frac{1}{8 + \frac{1}{2}}}}}$
Équation de Pell-Fermat
Soit $d$ un entier naturel qui n’est pas un carré. Le script suivant propose une méthode permettant de déterminer tous les couples $(x ; y)$ d’entiers naturels inférieurs à $n$ vérifiant l’équation de Pell-Fermat $x^2 − dy^2 = 1$.

  1. def pellfermat(d,n):
  2. result=[]
  3. for i in range(n+1):
  4. for j in range(1,n+1):
  5. if i**2==1+d*j**2:
  6. result.append((i,j))
  7. return result

Télécharger

Utilisation :

  1. >>> pellfermat(5,2000)
  2. [(9, 4), (161, 72)]

Télécharger

Et si l’on cherche celui dont la valeur de $y$ est minimale :

  1. def pellfermat2(d,n):
  2. result=[]
  3. for j in range(1,n+1):
  4. for i in range(1,n+1):
  5. if i**2==1+d*j**2:
  6. result.append((i,j))
  7. return result[0]

Télécharger

Utilisation :

  1. >>> pellfermat2(7,2000)
  2. (8, 3)

Télécharger

On peut alors généraliser et obtenir toutes les solutions de l’équation à l’aide de la suite $(x_n ; y_n), n \in \mathbb{N}^{{}*{}}$ définie par $x_n + y_n \sqrt{d} = \left(x_1 + y_1 \sqrt{d}\right)^n$ où $(x_1 ; y_1)$ est le couple de nombres obtenu par l’algorithme précédent, on a donc par exemple :
$(x_1 ; y_1) = (8 ; 3)$ et $x_2 + y_2 \sqrt{7} = \left(8 + 3 \sqrt{7}\right)^2 = 8^2 + 2\times 8\times 3\sqrt{7}+3^2\times 7 = 127 + 48\sqrt{7}$, donc $(x_2 ; y_2) = (127 ; 48)$.
Et en effet $127² - 7\times 48² = 1$.
Générateur de nombres aléatoires
Les nombres aléatoires sont presque toujours générés par des suites de nombres bien choisies (c’est-à-dire bien distribuées), voici le générateur de nombres aléatoires de votre calculatrice :

  1. suite = []
  2. def lehmer1(n):
  3. if n==0:
  4. return 12345
  5. else :
  6. return (lehmer1(n-1)*40014)%(2**31-85)
  7. def lehmer2(n):
  8. if n==0:
  9. return 67890
  10. else:
  11. return (lehmer2(n-1)*40692)%(2**31-249)
  12. def alea(n):
  13. alea= (lehmer1(n)-lehmer2(n))%(2**31-86)
  14. x=alea/(2**31-86)
  15. suite.append(x)
  16. return x

Télécharger

Utilisation :

  1. >>> alea(1)
  2. 0.9435974024931791

Télécharger

Pour vérifier que vous avez bien le même générateur dans votre calculatrice, saisissez sur celle-ci 0->NbrAléat (Sur TI, équivalent sur Casio), puis demandez un nombre aléatoire NbrAléat, vous obtenez :
0.9435974025
Soit le même nombre que ci-dessus (avec moins de chiffres sur la calculatrice), et vous pouvez continuer, les nombres suivants seront les mêmes (Voir l’article ici).
Pour avoir un aperçu graphique de la distribution aléatoire des termes de la suite :

  1. from tkinter import *
  2. def graphe(n):
  3. #n est le nombre de nombres aléatoires à générer
  4. fen = Tk()
  5. taille=200
  6. can = Canvas(fen, width =taille, height =taille, bg ='ivory')
  7. can.pack(side =TOP, padx =5, pady =5)
  8. for i in range(1,n):
  9. round(alea(i),10)
  10. #print(round(alea(i),10))
  11. zoom_x=taille/len(suite)
  12. zoom_y=taille/max(suite)
  13. for i in range(n-1):
  14. can.create_line(i*zoom_x-1, taille-zoom_y*suite[i]-1, i*zoom_x+1, taille-zoom_y*suite[i]+1, width=3, fill="blue")
  15. can.create_line(i*zoom_x-1, taille-zoom_y*suite[i]+1, i*zoom_x+1, taille-zoom_y*suite[i]-1, width=3, fill="blue")

Télécharger

Utilisation (Ne cherchez pas à demander trop de points ...) :

  1. >>> graphe(1000)


Changement de base
Nous allons nous intéresser dans ce paragraphe au changement d’ècriture d’un nombre en base a dans une base b.
Cela existe dans Python pour des bases particulières :
Pour convertir un nombre x (écrit sous forme d’une chaîne de caractères) d’une base quelconque en base 10 (ou nombre entier), on utilise la fonction int("x" ,base) où base représente la base sous laquelle est écrit le nombre :

  1. >>> int("AB", 16)
  2. 171
  3. >>> int("0xABB", 16)
  4. 2747
  5. >>> int("571", 8)
  6. 377
  7. >>> int("0571", 8)
  8. 377
  9. >>> int("1001010111010",2)
  10. 4794
  11. >>> int("0b100101011",2)
  12. 299

Télécharger

La méthode bin() convertit et renvoie l’équivalent binaire d’un nombre entier :

  1. >>> bin(23)
  2. '0b10111'

Télécharger

La fonction hex() convertit un nombre entier number en nombre hexadecimal qui lui correspond.

  1. >>> hex(255)
  2. '0xff'

Télécharger

Pour les autres conversions, l’on doit utiliser une fonction Python que l’on doit programmer, en voici une pour convertir un nombre entier décimal dans une base $p <10$ :

  1. def decimalversbase(n,b):
  2. res=""
  3. memoire=n
  4. while memoire>0:
  5. res=str(memoire%b)+res
  6. memoire=memoire//b
  7. return res

Télécharger

  1. >>> decimalversbase(1000,3)
  2. '1101001'
  3. >>> decimalversbase(1000,8)
  4. '1750'
  5. >>> decimalversbase(1000,2)
  6. '1111101000'

Télécharger

Et une autre pour convertir un nombre dans une base $p <10$ en un entier décimal :

  1. def baseversdecimal(n,b):
  2. res=0
  3. if n<b:
  4. res=n
  5. else:
  6. n=str(n)
  7. l=len(str(n))
  8. for i in range(l):
  9. res+=int(n[l-i-1])*b**i
  10. return res

Télécharger

Utilisation :

  1. >>> baseversdecimal(1101001,3)
  2. 1000
  3. >>> baseversdecimal(1750,8)
  4. 1000

Télécharger

Et en combinant les deux on obtient la conversion d’un nombre en base $a < 10$ en un nombre en base $b < 10$ :

  1. def baseaversbaseb(n,a,b):
  2. return decimalversbase(baseversdecimal(n,a),b)

Télécharger

Utilisation :

  1. >>> baseaversbaseb(1750,8,3)
  2. '1101001'

Télécharger

Pour des conversions dans des bases inférieures à 36 :

  1. from math import pow, log
  2. def convert(n,a,b):
  3. n=str(n)
  4. " n est un str, a et b des entiers "
  5. chiffres = ['0','1','2','3','4','5','6','7','8','9','a','b','c','d','e','f','g','h','i','j','k','l','m','n','o','p','q','r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z']
  6. nb_decimal = 0
  7. i = len(n)-1
  8. for x in n:
  9. nb_decimal += chiffres.index(x)*int(pow(a,i))
  10. i= i-1
  11. nbr_chars = int(log(nb_decimal)/log(b)+1)
  12. res = ""
  13. for i in range(0, nbr_chars):
  14. x = int(nb_decimal / pow(b,nbr_chars-1-i))
  15. nb_decimal -= int(x * pow(b,nbr_chars-1-i))
  16. res = res+chiffres[x]
  17. return res

Télécharger

Utilisation :

  1. >>> convert(1000, 12, 5)
  2. '23403'
  3. >>> convert(1000, 35, 5)
  4. '2333000'

Télécharger

Graphes et matrices

MATRICES
Calculs de base
Nous sommes tentés de faire cela :

  1. >>> A=[[2,-1,3],[4,-2,0],[-3,3,6]]
  2. >>> B=[[ 1,3,11],[7,9,17],[13,15,5]]
  3. >>> 3*A-2*B
  4. Traceback (most recent call last):
  5. File "<pyshell#37>", line 1, in <module>
  6. 3*A-2*B
  7. TypeError: unsupported operand type(s) for -: 'list' and 'list'

Télécharger

Cela ne fonctionne pas parce que A et B ne sont pas des matrices pour Python, rappelez-vous que l’on a plusieurs fois utilisé cela :

  1. >>> A=[0]
  2. >>> 10*A
  3. [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

Télécharger

Forcément le calcul de 3*A-2*B allait poser problème ...
Une solution est d’utiliser le module numpy :

  1. >>> import numpy as np
  2. >>> A=[[2,-1,3],[4,-2,0],[-3,3,6]]
  3. >>> A
  4. [[2, -1, 3], [4, -2, 0], [-3, 3, 6]]
  5. >>> A=np.array(A)
  6. >>> A
  7. array([[ 2, -1, 3],
  8. [ 4, -2, 0],
  9. [-3, 3, 6]])

Télécharger

La commande array() tranforme la liste de listes A en un tableau A.
On peut maintenant calculer :

  1. >>> B=np.array([[ 1,3,11],[7,9,17],[13,15,5]])
  2. >>> 3*A-2*B
  3. array([[ 4, -9, -13],
  4. [ -2, -24, -34],
  5. [-35, -21, 8]])
  6. >>> -A
  7. array([[-2, 1, -3],
  8. [-4, 2, 0],
  9. [ 3, -3, -6]])

Télécharger

Mais cela s’arrète là ! Impossible de calculer 1/A ou A³, enfin si mais les calculs ne sont pas des calculs matriciels :

  1. >>> 1/A
  2. Warning (from warnings module):
  3. File "<pyshell#50>", line 1
  4. RuntimeWarning: divide by zero encountered in true_divide
  5. array([[ 0.5 , -1. , 0.33333333],
  6. [ 0.25 , -0.5 , inf],
  7. [-0.33333333, 0.33333333, 0.16666667]])
  8. >>> A**2
  9. array([[ 4, 1, 9],
  10. [16, 4, 0],
  11. [ 9, 9, 36]])

Télécharger

En fait 1/A renvoie le tableau des coefficients de A inversés et A**2 renvoie le tableau des coefficients de A au carré ...
numpy permet cependant ces calculs en utilisant des commandes du module, pour l’inverse :

  1. >>> np.linalg.inv(A)
  2. array([[-0.66666667, 0.83333333, 0.33333333],
  3. [-1.33333333, 1.16666667, 0.66666667],
  4. [ 0.33333333, -0.16666667, 0. ]])

Télécharger

Si l’on veut calculer le produit matriciel A.B :

  1. >>> np.dot(A,B)
  2. array([[ 34, 42, 20],
  3. [-10, -6, 10],
  4. [ 96, 108, 48]])

Télécharger

Autre possibilité :

  1. >>> A.dot(B)
  2. array([[ 34, 42, 20],
  3. [-10, -6, 10],
  4. [ 96, 108, 48]])

Télécharger

Si l’on veut calculer le cube de A :

  1. >>> np.linalg.matrix_power(A,3)
  2. array([[-54, 63, 117],
  3. [-36, 36, 72],
  4. [-45, 63, 126]])

Télécharger

Solutions d’un système d’équations linéaires
L’opération peut-être la plus fréquente en algèbre linéaire est la solution de systèmes d’équations algébriques linéaires : $Ax = b$, où A est la matrice des coefficients, b est un vecteur donné et x est le vecteur solution. La fonction np.linalg.solve (A, b) fait le travail :
Pour résoudre ce système :
$ \left\{ \begin{array}{r c l} 2 x + y &=& 15\\ -x + 4 y &=&-7 \end{array} \right. $

  1. >>> A = np.array([[2, 1], [-1, 4]])
  2. >>> b= np.array([15, -7])
  3. >>> np.linalg.solve(A, b)
  4. array([7.44444444, 0.11111111])

Télécharger

On peut bien sûr se passer de cette commande puisque si A est inversible $Ax = b \Leftrightarrow x = A^{-1}b$ :

  1. >>> np.dot(np.linalg.inv(A),b)
  2. array([7.44444444, 0.11111111])

Télécharger

Inversion d’une matrice M carrée d’ordre 2
On souhaite inverser la matrice $M = \begin{pmatrix} m_{11} & m_{12}\\ m_{21} & m_{22} \end{pmatrix}$

  1. def inversion2(m11,m12,m21,m22):
  2. d=m11*m22-m12*m21
  3. if d==0:
  4. return "La matrice M n'est pas inversible car son déterminant est nul."
  5. else:
  6. return [[m22/d,-m12/d],[-m21/d,m11/d]]

Télécharger

Utilisation :

  1. >>> inversion2(1,0,-2,3)
  2. [[1.0, 0.0], [0.6666666666666666, 0.3333333333333333]]
  3. >>> inversion2(1,-1,-2,3)
  4. [[3.0, 1.0], [2.0, 1.0]]
  5. >>> inversion2(1,-1,-3,3)
  6. "La matrice M n'est pas inversible car son déterminant est nul."

Télécharger

Matrices particulières :

Commande Description (rapide)
np.zeros(n) ((n,p)) vecteur nul de taille n, matrice nulle de taille (n,p)
np.eye(n) (n,p) matrice de taille n ou (n,p) avec des 1 sur la diagonale et des zéros ailleurs
np.ones(n) ((n,p)) vecteur de taille n, matrice de taille n,p remplie de 1
np.diag(v) matrice diagonale dont la diagonale est le vecteur v
np.diag(v,k)         matrice dont la ‘diagonale’ décalée de k est le vecteur v (k est un entier relatif)
np.random.rand(n) (n,p)         vecteur (taille n), matrice (taille n,p) à coefficients aléatoires uniformes sur [0,1]

Utilisation :

  1. >>> np.zeros(3)
  2. array([0., 0., 0.])
  3. >>> np.zeros((3,3))
  4. array([[0., 0., 0.],
  5. [0., 0., 0.],
  6. [0., 0., 0.]])
  7. >>> np.zeros((3,2))
  8. array([[0., 0.],
  9. [0., 0.],
  10. [0., 0.]])
  11. >>> np.eye(3)
  12. array([[1., 0., 0.],
  13. [0., 1., 0.],
  14. [0., 0., 1.]])
  15. >>> np.eye(3,2)
  16. array([[1., 0.],
  17. [0., 1.],
  18. [0., 0.]])
  19. >>> np.ones(3)
  20. array([1., 1., 1.])
  21. >>> np.ones((3,2))
  22. array([[1., 1.],
  23. [1., 1.],
  24. [1., 1.]])
  25. >>> np.diag([-2,3,1])
  26. array([[-2, 0, 0],
  27. [ 0, 3, 0],
  28. [ 0, 0, 1]])
  29. >>> np.diag([-2,3,1],1)
  30. array([[ 0, -2, 0, 0],
  31. [ 0, 0, 3, 0],
  32. [ 0, 0, 0, 1],
  33. [ 0, 0, 0, 0]])
  34. >>> np.diag([-2,3,1],-2)
  35. array([[ 0, 0, 0, 0, 0],
  36. [ 0, 0, 0, 0, 0],
  37. [-2, 0, 0, 0, 0],
  38. [ 0, 3, 0, 0, 0],
  39. [ 0, 0, 1, 0, 0]])
  40. >>> np.random.rand(5)
  41. array([0.5704163 , 0.46493989, 0.85085904, 0.56177765, 0.06711559])
  42. >>> np.random.rand(5,2)
  43. array([[0.69873188, 0.66611399],
  44. [0.76550507, 0.25193332],
  45. [0.28569524, 0.67081913],
  46. [0.63928606, 0.97949063],
  47. [0.30312998, 0.31526063]])

Télécharger

GRAPHES
Distribution invariante d’une chaîne de Markov

On souhaite déterminer la distribution invariante $\pi(a,b)$ d’une chaîne de Markov à deux états dont la matrice de transition est $ P = \begin{pmatrix} 1 - \lambda & \lambda\\ \mu & 1 - \mu \end{pmatrix} $. Or $\pi P = \pi \Leftrightarrow \left\{ \begin{array}{r c l} -\lambda a + b\mu &=& 0\\ a + b &=& 1 \end{array} \right. $, on a donc une équation $A\pi = B$ et donc $\pi = A^{-1}B$ d’où l’algorithme :

  1. import numpy as np
  2. def distribution_invariante(lambd,mu):
  3. A=np.array([[-lambd,mu],[1,1]])
  4. B=np.array([0,1])
  5. pi=np.dot(np.linalg.inv(A),B)
  6. return pi

Télécharger

Utilisation :

  1. >>> distribution_invariante(0.4,0.7)
  2. array([0.63636364, 0.36363636])

Télécharger

Vérification (On regarde si $\pi P = \pi$) :

  1. >>> np.dot([0.63636364, 0.36363636],[[0.6,0.4],[0.7,0.3]])
  2. array([0.63636364, 0.36363636])

Télécharger

Le modèle des urnes d’Ehrenfest
On considère deux urnes 1 et 2, ainsi que k boules, numérotées de 1 à k. Initialement, toutes les boules se trouvent dans l’urne 1. Le processus consiste à répéter l’opération suivante :
- Tirer au hasard un numéro i compris entre 1 et k, prendre la boule n°i, la transférer dans l’urne où elle n’était pas.
Dans Python, les boules sont numérotées de 0 à k-1, donc il y en a toujours k :

  1. from random import randint
  2. def ehrenfest(n,k):
  3. L=[1]*k #On met les k boules dans l'urne 1
  4. for i in range(n):
  5. choix=randint(0,k-1) #On choisit une boule au hasard
  6. if L[choix]==1:
  7. L[choix]=2 #Si elle est dans l'urne 1, on la met dans la 2
  8. else:
  9. L[choix]=1 #Si elle est dans l'urne 2, on la met dans la 1
  10. return L

Télécharger

Utilisation :

  1. >>> ehrenfest(1000,10)
  2. [1, 1, 1, 1, 2, 1, 2, 2, 1, 2]
  3. >>> ehrenfest(1000,15)
  4. [2, 2, 1, 1, 1, 2, 1, 2, 1, 1, 2, 1, 1, 2, 1]

Télécharger

Dans ce modèle, on suit au cours du temps t (discret) le nombre total de boules n(t) présentes dans l’urne 1. On obtient une courbe qui part initialement de n(0)=k et commence par décroître vers la valeur moyenne k/2.

  1. def graphique(n,k):
  2. L=[1]*k
  3. x=[0]
  4. y=[k]
  5. z=[k//2]*(n+1)
  6. nb=k
  7. for i in range(n):
  8. choix=randint(0,k-1) #On choisit une boule au hasard
  9. x.append(i)
  10. if L[choix]==1:
  11. L[choix]=2 #Si elle est dans l'urne 1, on la met dans la 2
  12. nb=nb-1
  13. else:
  14. L[choix]=1 #Si elle est dans l'urne 2, on la met dans la 1
  15. nb=nb+1
  16. y.append(nb)
  17.  
  18. pl.plot(x,y)
  19. pl.plot(x,z)
  20. pl.show()
  21. return L.count(1)/k

Télécharger

Utilisation :

  1. >>> graphique(100,4)
  2. 0.5

Télécharger

k = 4 boules ; 100 tirages. Les fluctuations autour de la moyenne sont importantes lorsque k est petit, et les récurrences à l’état initial sont particulièrement visibles.

Mais cette décroissance est irrégulière : il existe des fluctuations autour de la valeur moyenne k/2, qui peuvent devenir parfois très importantes (ceci est particulièrement visible lorsque k est petit).
En particulier, quel que soit le nombre de boules k fini, il existe toujours des récurrences à l’état initial, pour lesquelles toutes les boules reviennent dans l’urne 1 après une durée finie. Mais, comme le temps moyen entre deux récurrences consécutives croît très rapidement avec k, ces récurrences ne nous apparaissent pas lorsque k est très grand.

  1. >>> graphique(10000,500)
  2. 0.52

Télécharger


Construction de graphes
Avec le module Graphviz, pour obtenir le graphe orienté de la chaine de Markov à deux états vu précédemment :

  1. from graphviz import Digraph
  2. g = Digraph('G', filename='Markov', engine='dot')
  3. g.edge("0", "0",label="1 - μ")
  4. g.edge("1", "0",label="λ")
  5. g.edge("1", "1",label="1 - λ")
  6. g.edge("0", "1",label="μ")
  7. g.format="svg"
  8. g.render(view=True)

Télécharger

Utilisation :

  1. from graphviz import Digraph
  2. g = Digraph('G', filename='Distrib', engine='dot')
  3. g.edge("E", "E",label="0.3")
  4. g.edge("A", "E",label="0.4")
  5. g.edge("A", "A",label="0.6")
  6. g.edge("E", "A",label="0.7")
  7. g.format="svg"
  8. g.render(view=True)

Télécharger

Utilisation :

Pour un graphe non orienté, pondéré :

  1. from graphviz import Graph
  2. g = Graph('G', filename='Itineraire', engine='circo')
  3. g.edge("Montpellier", "Béziers",label="72 km, 57 min")
  4. g.edge("Montpellier", "Nîmes",label="56 km, 50 min")
  5. g.edge("Alès", "Nîmes",label="55 km, 38 min")
  6. g.edge("Ganges", "Alès",label="51 km, 55 min")
  7. g.edge("Ganges", "Nîmes",label="67 km, 68 min")
  8. g.edge("Ganges", "Béziers",label="98 km, 89 min")
  9. g.edge("Ganges", "Montpellier",label="47 km, 54 min")
  10. g.edge("Alès", "Montpellier",label="80 km, 96 min")
  11. g.format="svg"
  12. g.render(view=True)

Télécharger

Utilisation :

On peut obtenir un rendu différent avec le paramètre engine (Voir paragraphe Roadmap en suivant ce lien) :

  1. g = Graph('G', filename='Itineraire', engine='dot')
  1. g = Graph('G', filename='Itineraire', engine='sfdp')

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