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 :
- def pgcd(a,b):
- if b==0:
- return a
- else:
- while a!=b:
- if a>=b:
- a -= b
- else:
- b -= a
- return a
Puis l’algorithme des divisions successives :
- def pgcd(a,b):
- if b==0:
- return a
- else:
- while b!=0:
- a,b=b,a%b
- return a
Le même en version récursive :
- def pgcd(a,b):
- if b==0:
- return a
- else:
- return pgcd(b,a%b)
Ce même algorithme peut être utilisé pour tester si deux nombres sont premiers entre eux :
- def premier(a,b):
- while b!=0:
- r=a%b
- a=b
- b=r
- return a==1
Calcul d’un couple de Bézout :
- def bezout(a,b):
- (u0,v0,u1,v1)=(1,0,0,1)
- while b:
- (q,r)=divmod(a,b)
- (a,b)=(b,r)
- (u0,v0,u1,v1)=(u1,v1,u0-q*u1,v0-q*v1)
- return (u0,v0)
Utilisation :
- >>> bezout(412,326)
- (-72, 91)
Version récursive :
- def bezout(a, b):
- if b == 0:
- return (1,0)
- else:
- (u,v) = bezout(b,a%b)
- return (v, u-(a//b)*v)
On peut aussi calculer ce couple en itérant sur u tant que v =(1 - au)/b n’est pas entier :
- def Bezout(a,b):
- if premier(a,b)==False:
- return "Pas de solution ! Les nombres ne sont pas premiers entre eux !"
- else:
- u=0
- v=1/b
- while v!=int(v):
- u+=1
- v=(1-a*u)/b
- return u,int(v)
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 :
- def eratosthene(n):
- """retourne le tableau des nombres premiers <= n (crible d'eratosthene)"""
- n += 1
- tableau = [0,0]+[i for i in range(2,n)]
- for i in range(2,n):
- if tableau[i]!= 0:
- #c'est un nombre 1er : on le garde, mais on supprime ses multiples
- for j in range(i*2,n,i):
- tableau[j] = 0
- return [p for p in tableau if p!=0]
On peut améliorer ce code (je n’en suis plus très sûr maintenant mais je laisse ça là ...) :
- def eratosthene2(n):
- """retourne la liste des nombres premiers <= n (crible d'eratosthene)"""
- if n<2:
- return []
- n += 1
- tableau = [False,False] + [True]*(n-2)
- tableau[2::2] = [False]*((n-2)//2 + n%2) # sup. des nb pairs
- premiers = [2] # initialisation du tableau des nb 1ers (2 est 1er)
- racine = int(n**0.5)
- racine = racine + [1,0][racine%2] # pour que racine soit impair
- for i in range(3, racine+1, 2):
- if tableau[i]:
- # on enregistre le nouveau nb 1er
- premiers.append(i)
- # on élimine i et ses multiples
- tableau[i::i] = [False]*((n-i)//i + int((n-i)%i>0))
- for i in range(racine, n, 2):
- if tableau[i]:
- # on enregistre le nouveau nb 1er (pas de multiples à supprimer)
- premiers.append(i)
- return premiers
Ou en utilisant les listes en compréhension :
- def crible(N):
- crible=[i for i in range(2,N+1)]
- for n in range(2,round(sqrt(N))+1):
- if n in crible: #ce test sert à ne pas éliminer les multiples de 4 alors que l'on a déjà éliminé les multiples de 2.
- crible=[k for k in crible if k<=n or k%n!=0]
- return crible
Utilisation :
- >>> crible(100)
- [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]
Décomposition en facteurs (...)
Décomposition en facteurs premiers
- def facteurs(n):
- F=[]
- x=n
- i=2
- while i<n//2+1:
- a,b=divmod(x,i)
- if b==0:
- F.append(i)
- x=a
- else:
- i+=1
- if F==[]:
- F.append(n)
- return F
On peut accélérer les calculs en arrangeant un peu le code :
- from math import sqrt
- def facteurs(n):
- """facteurs(n): décomposition d'un nombre entier n en facteurs premiers"""
- F = []
- if n==1:
- return F
- # recherche de tous les facteurs 2 s'il y en a
- while n>=2:
- x,r = divmod(n,2)
- if r!=0:
- break
- F.append(2)
- n = x
- # recherche des facteurs 1er >2
- i=3
- rn = int(sqrt(n))
- while i<=n:
- if i>rn:
- F.append(n)
- break
- x,r = divmod(n,i)
- if r==0:
- F.append(i)
- n=x
- rn = int(sqrt(n))
- else:
- i += 2
- return F
Bernard Parisse, lecteur attentif de Mathematice et auteur de XCAS, nous propose ce script qui n’utilise pas divmod et que nous publions volontiers [2] :
- def facto(n):
- p=2
- l=[]
- while p*p<=n:
- if n%p==0:
- l.append(p)
- n=n//p
- else:
- p+=1
- if n!=1:
- l.append(n)
- return l
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 :
- from random import *
- def PR(n):
- nAccueil,n1,n2=0,0,0
- start=randint(0,2)
- if start==0:
- page="Accueil"
- nAccueil=1
- elif start==1:
- page="page 1"
- n1=1
- else:
- page="page 2"
- n2=1
- for i in range(1,n):
- if page=="Accueil":
- alea=randint(0,1)
- if alea==0:
- page="page 1"
- n1=n1+1
- else:
- page="page 2"
- n2=n2+1
- elif page=="page 1":
- page="Accueil"
- nAccueil=nAccueil+1
- else:
- alea=randint(0,1)
- if alea==0:
- page="Accueil"
- nAccueil=nAccueil+1
- else:
- page="page 1"
- n1=n1+1
- fAccueil,f1,f2=round(nAccueil/n,2),round(n1/n,2),round(n2/n,2)
- return fAccueil,f1,f2
Utilisation :
- >>> PR(100)
- (0.46, 0.33, 0.21)
- >>> PR(1000)
- (0.45, 0.34, 0.21)
- >>> PR(100000)
- (0.44, 0.33, 0.22)
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 ...) :
- >>> z1=complex(-2,3)
- >>> z1
- (-2+3j)
- >>> z1.real
- -2.0
- >>> z1.imag
- 3.0
- >>> z1.conjugate()
- (-2-3j)
- >>> z2=complex(1,-1)
- >>> z1+z2
- (-1+2j)
- >>> z1-z2
- (-3+4j)
- >>> z1*z2
- (1+5j)
- >>> z1/z2
- (-2.5+0.5j)
- >>> z2**4
- (-4-0j)
- >>> 1/complex(1,1)
- (0.5-0.5j)
- >>> abs(complex(-1,0))
- 1.0
- >>> abs(complex(1,1))
- 1.4142135623730951
abs(z)
renvoie le module du nombre complexe z.
Et en utilisant le module cmath
(Fonctions mathématiques pour nombres complexes) :
- from cmath import *
- >>> phase(z1)
- 2.158798930342464
- >>> phase(complex(-1,0))
- 3.141592653589793
- >>> phase(complex(0,1))
- 1.5707963267948966
- >>> polar(complex(10,0))
- (10.0, 0.0)
- >>> polar(complex(0,-1))
- (1.0, -1.5707963267948966)
- >>> polar(complex(-1,0))
- (1.0, 3.141592653589793)
- >>> rect(2,pi/2)
- (1.2246467991473532e-16+2j)
- >>> rect(2,pi/3)
- (1.0000000000000002+1.7320508075688772j)
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]$ :
- def arg(z):
- r=abs(z)
- c=z.real/r
- if z.imag>=0:
- a=math.acos(c)
- else:
- a=-math.acos(c)
- return a
Utilisation et vérification :
- >>> arg(complex(1,-1))
- -0.7853981633974484
- >>> arg(complex(1,1))
- 0.7853981633974484
- >>> arg(complex(-1,1))
- 2.356194490192345
- >>> arg(complex(-1,-1))
- -2.356194490192345
- >>> phase(complex(1,-1))
- -0.7853981633974483
- >>> phase(complex(1,1))
- 0.7853981633974483
- >>> phase(complex(-1,-1))
- -2.356194490192345
- >>> phase(complex(-1,1))
- 2.356194490192345
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 :
- def sqrt_complex(c):
- r=sqrt(abs(c))
- phi=phase(c)/2
- x=r*cos(phi)
- y=r*sin(phi)
- return complex(-x,-y),complex(x,y)
Utilisation :
- >>> sqrt_complex(complex(0,1))
- ((-0.7071067811865476-0.7071067811865475j), (0.7071067811865476+0.7071067811865475j))
- >>> sqrt_complex(complex(1,1))
- ((-1.0986841134678098-0.45508986056222733j), (1.0986841134678098+0.45508986056222733j))
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 :
- def sqrt_complex2(c):
- a=c.real
- b=c.imag
- module=abs(c)
- x=sqrt((a+module)/2)
- y=sqrt((module-a)/2)
- if b<0:
- y=-y
- return complex(-x,-y),complex(x,y)
Solutions dans C d’une équation du second degré à coefficients réels
- def polynome_degre2(a,b,c):
- delta=b**2-4*a*c
- if delta>0:
- print("1")
- return ["Deux solutions réelles :",(-b-math.sqrt(delta))/(2*a),(-b+math.sqrt(delta))/(2*a)]
- elif b==0:
- print("2")
- return ["Une solution réelle :",-b/(2*a)]
- else:
- print("3")
- return ["Deux solutions complexes :",complex(-b/(2*a),-math.sqrt(-delta)/(2*a)),complex(-b/(2*a),math.sqrt(-delta)/(2*a))]
Utilisation :
- >>> polynome_degre2(1,1,0)
- ['Deux solutions réelles :', -1.0, 0.0]
- >>> polynome_degre2(25,-4,1/4)
- ['Deux solutions complexes :', (0.08-0.06j), (0.08+0.06j)]
- >>> polynome_degre2(1,-14,49)
- ['Une solution réelle :', 7.0]
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.
- def Horner(L,a):
- n=len(L)
- b=L[0]
- for i in range(1,n):
- b=b*a+L[i]
- return b
Utilisation (avec des complexes ou des réels) :
- >>> Horner([2,-5,1,3],10)
- 1513
- >>> Horner([2,-5,1,3],1j)
- (8-1j)
- >>> Horner([2,-5,1,3],complex(1,1))
- -5j
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 :
- from pylab import *
- from cmath import *
- def polygone(n):
- for k in range(n):
- z=2*exp(1j*2*k*pi/n)
- plot(z.real,z.imag,"r.",markersize=12)
- show()
Utilisation :
- >>> 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 :
- from cmath import *
- def Points_alignes(a,b,c):
- z=(c-a)/(b-a)
- if phase(z)==0 or phase(z)==pi:
- return True
- else:
- return False
Utilisation :
- >>> Points_alignes(3+5j,-4+4j,-11+3j)
- True
- >>> Points_alignes(1j,8,-7+2j)
- False
- >>> Points_alignes(1-1j,5+1j,-1-2j)
- True
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 :
- from cmath import *
- def Droites_orthogonales(a,b,c):
- z=(c-a)/(b-a)
- if phase(z)==pi/2 or phase(z)==-pi/2:
- return True
- else:
- return False
>>> 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.
- import matplotlib.pyplot as plt
- from random import random
- def mandelbrot(n):
- for i in range(n):
- c=complex(4*random()-2,4*random()-2)
- z=0
- k=0
- while abs(z)<2 and k<50:
- z=z**2+c
- k=k+1
- if k==50:
- plt.plot(c.real,c.imag,"rx",ms=0.5)
- plt.xlim=(-2,2)
- plt.ylim=(-2,2)
- ax=plt.gca()
- ax.spines["bottom"].set_position(("data",0))
- ax.spines["left"].set_position(("data",0))
- plt.show()
Utilisation :
- >>> mandelbrot(100000)
- >>> mandelbrot(500000)
- >>> 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$.
- from PIL import Image
- def Julia(c,taille):
- a=c.real
- b=c.imag
- xmin = -1.25
- xmax = 1.25
- ymin = -1.25
- ymax = 1.25
- iterationmax = 200
- im = Image.new("RGB",(taille,taille),(255,255,255))
- pixels = im.load()
- for line in range(taille):
- for col in range(taille):
- i = 1
- x = xmin+col*(xmax-xmin)/taille
- y = ymax-line*(ymax-ymin)/taille
- while i<=iterationmax and (x**2+y**2)<=4:
- stock = x
- x = x**2-y**2+a
- y = 2*stock*y+b
- i += 1
- if i>iterationmax and (x**2+y**2)<=4:
- pixels[col,line] = (0,0,0)
- else:
- pixels[col,line] = ((4*i)%256,2*i,(6*i)%256)
- im.save("julia1.png")
- im.show()
Utilisation (Pensez à changer le nom du fichier sauvegardé si vous ne voulez pas perdre l’image précédente) :
- >>> Julia(complex(0.284, 0.0122),400)
- >>> Julia(complex(-0.1011,0.9563),400)
- Julia(complex(-0.1528,1.0397),400)
- Julia(complex(-0.06783611264225832,0.6617460391250546),400)
- 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.
- >>> Julia(complex(-0.12256,0.74486),400)
Puissance n-ième d’un nombre complexe
On peut se servir du calcul des coefficients binomiaux vu en Spécialité pour calculer la puissance n-ième d’un nombre complexe en utilisant le binôme de Newton. binomialp(n,p) calcule $\left(\begin{matrix}
n\\
p
\end{matrix}\right)$ en utilisant la relation de récurrence $\left(\begin{matrix}
n\\
p
\end{matrix}\right) = \left(\begin{matrix}
n-1\\
p-1
\end{matrix}\right) + \left(\begin{matrix}
n-1\\
p
\end{matrix}\right) $ et binomial(n) renvoie la liste des coefficients binomiaux au rang n.
- def binome(n):
- pascal=[1,0]
- k=1
- while k<=n:
- k=k+1
- pascal2=[1]
- for i in range(1,k):
- pascal2.append(pascal[i-1]+pascal[i])
- pascal=pascal2
- pascal.append(0)
- pascal.remove(0)
- return pascal
- def developpe(a,b,k):
- S=0
- L=binome(k)
- for i in range(k+1):
- S=S+L[i]*a**(k-i)*b**i
- return S
Utilisation :
- >>> developpe(5,complex(0,3),5)
- (-6100+2868j)
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.
- def equation_diophantienne(a,b,c):
- p=pgcd(a,b)
- x0=c*bezout(a,b)[0]/p
- y0=c*bezout(a,b)[1]/p
- return [(x0,y0),(b/p,"k +",x0,";"-a/p,"k +",y0)]
Utilisation
- >>> equation_diophantienne(3,2,12)
- [(12.0, -12.0), (2.0, 'k +', 12.0, ';', -3.0, 'k +', -12.0)]
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 :
- def inverse_modulaire(a,m):
- if pgcd(a,m)>1:
- return None
- else:
- return m+bezout(a,m)[0]
Utilisation :
- >>> inverse_modulaire(6,31)
- 26
- >>> inverse_modulaire(6,7)
- 6
- >>> inverse_modulaire(4,8)
- >>> inverse_modulaire(13,17)
- 21
- >>> inverse_modulaire(39,53)
- 34
On peut aussi déterminer cet inverse modulaire sans Bézout :
- def inverse_modulaire(x,m):
- for n in range(m):
- if (x*n)%m == 1:
- return n
- elif n == m-1:
- return "Null"
- else:
- continue
Diviseurs
Si l’on veut afficher la liste des diviseurs d’un nombre dans $\mathbb{N}$ :
- def diviseurs(n):
- L=[]
- for i in range(1,n+1):
- if n%i==0:
- L.append(i)
- return L
Utilisation :
- >>> diviseurs(4567)
- [1, 4567]
- >>> diviseurs(4568)
- [1, 2, 4, 8, 571, 1142, 2284, 4568]
- >>> diviseurs(60)
- [1, 2, 3, 4, 5, 6, 10, 12, 15, 20, 30, 60]
- >>> diviseurs(36)
- [1, 2, 3, 4, 6, 9, 12, 18, 36]
Avec une liste en compréhension :
- def div(n):
- return [x for x in range(1,n+1) if n%x==0]
Utilisation :
- >>> div(123456789)
- [1, 3, 9, 3607, 3803, 10821, 11409, 32463, 34227, 13717421, 41152263, 123456789]
Si l’on veut afficher la liste des diviseurs d’un nombre dans $\mathbb{Z}$ :
- def diviseursZ(n):
- if n<0:
- M=diviseurs(-n)
- else:
- M=diviseurs(n)
- L=[]
- for x in M:
- L.append(-x)
- L=L+M
- L.sort()
- return L
Utilisation :
- >>> diviseursZ(72)
- [-72, -36, -24, -18, -12, -9, -8, -6, -4, -3, -2, -1, 1, 2, 3, 4, 6, 8, 9, 12, 18, 24, 36, 72]
- >>> diviseursZ(-72)
- [-72, -36, -24, -18, -12, -9, -8, -6, -4, -3, -2, -1, 1, 2, 3, 4, 6, 8, 9, 12, 18, 24, 36, 72]
Avec une liste en compréhension :
- def divZ(n):
- return [i for i in [x for x in range(-n,n+1) if x!=0] if n%i==0]
Utilisation :
- >>> divZ(124)
- [-124, -62, -31, -4, -2, -1, 1, 2, 4, 31, 62, 124]
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 :
- def est_premier(n):
- L=diviseurs(n)
- if len(L)==2:
- return True
- else:
- return False
Ou plus simplement :
- def est_premier(n):
- L=diviseurs(n)
- return len(L)==2
Utilisation :
- >>> est_premier(1)
- False
- >>> est_premier(2)
- True
- >>> est_premier(0)
- False
- >>> est_premier(6)
- False
- >>> est_premier(19)
- True
Liste des n premiers nombres premiers
- def liste_premiers(n):
- i=0
- j=2
- premiers=[]
- while i<n:
- if est_premier(j):
- premiers.append(j)
- j=j+1
- i=i+1
- else:
- j=j+1
- return premiers
Utilisation :
- >>> liste_premiers(25)
- [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]
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 :
- import matplotlib.pyplot as plt
- def graphe(m,n):
- for colonne in range(1,m+1):
- for ligne in range(1,n+1):
- nombre =colonne+(ligne-1)*m
- if est_premier(nombre)==True:
- plt.plot(ligne,colonne,"ko")
- plt.xticks([],[])
- plt.yticks([],[])
- plt.show()
Utilisation :
- >>> 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.")
:
- >>> graphe(100,100)
Nombres amicaux
On appelle diviseur strict de $n$ un nombre entier naturel qui divise $n$, strictement inférieur à $n$ :
- def diviseurs_stricts(n):
- L=[]
- for i in range(1,n//2+1):
- if n%i==0:
- L.append(i)
- return L
Utilisation :
- >>> diviseurs_stricts(220)
- [1, 2, 4, 5, 10, 11, 20, 22, 44, 55, 110]
- >>> diviseurs_stricts(284)
- [1, 2, 4, 71, 142]
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 :
- def amicaux(n,p):
- if sum(diviseurs_stricts(n))==p and sum(diviseurs_stricts(p))==n:
- return True
- else:
- return False
Utilisation :
- >>> amicaux(36,48)
- False
- >>> amicaux(220,284)
- True
On peut alors chercher tous les amis inférieurs à n :
- def cherche_amis(n):
- amis=[]
- for i in range(2,n):
- for j in range(i+1,n):
- if amicaux(i,j)==True:
- amis.append((i,j))
- return amis
Et ils ne sont pas nombreux :
- >>> cherche_amis(1000)
- [(220, 284)]
- >>> cherche_amis(2000)
- [(220, 284), (1184, 1210)]
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 :
- def est_parfait(n):
- if sum(diviseurs_stricts(n))==n:
- return True
- else:
- return False
Utilisation :
- >>> est_parfait(6)
- True
- >>> est_parfait(28)
- True
- >>> est_parfait(36)
- False
On peut alors chercher tous les parfaits inférieurs à n :
- def cherche_parfait(n):
- parfaits=[]
- for i in range(2,n):
- if est_parfait(i)==True:
- parfaits.append(i)
- return parfaits
Utilisation :
- >>> cherche_parfait(500)
- [6, 28, 496]
- >>> cherche_parfait(5000)
- [6, 28, 496]
- >>> cherche_parfait(50000)
- [6, 28, 496, 8128]
On peut faire cette recherche directement, sans utiliser de fonctions annexes :
- def parfaits(n):
- L=[]
- for i in range(2,n+1):
- s=0
- for j in range(1,i//2+1):
- if i%j==0:
- s+=j
- if s==i:
- L.append(i)
- return L
Utilisation :
- >>> parfaits(30)
- [6, 24, 28]
- >>> parfaits(500)
- [6, 24, 28, 496]
Tableau de congruences
En se servant du script ci-dessous, l’on peut, sans trop de difficulté obtenir un tableau de congruences :
- print("{:22}".format("n ≡ ... [4] |"),end="")
- for n in range(4):
- print("{:2d}".format(n%4),end=" | ")
- print()
- print("{:22}".format("2n ≡ ... [4] | "),end="")
- for n in range(4):
- print("{:2d}".format(2*n %4),end=" | ")
- print()
- print("{:22}".format("n² + 1 ≡ ... [4] | "),end="")
- for n in range(4):
- print("{:2d}".format((n**2+1)%4),end=" | ")
- print()
- print("{:22}".format("2n(n² + 1) ≡ ... [4]| "),end="")
- for n in range(4):
- print("{:2d}".format((2*n*(n**2+1))%4),end=" | ")
Utilisation :
- n ≡ ... [4] | 0 | 1 | 2 | 3 |
- 2n ≡ ... [4] | 0 | 2 | 0 | 2 |
- n² + 1 ≡ ... [4] | 1 | 2 | 1 | 2 |
- 2n(n² + 1) ≡ ... [4]| 0 | 0 | 0 | 0 |
Décomposition de la fraction $\frac{n}{d}$ en fraction continue
- def frac_continue(n,d):
- solution=[]
- while d:
- a=n//d
- solution.append(a)
- n,d=d,n-a*d
- return solution
Utilisation :
- >>> frac_continue(314,100)
- [3, 7, 7]
- >>> frac_continue(31415,10000)
- [3, 7, 14, 1, 8, 2]
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$.
- def pellfermat(d,n):
- result=[]
- for i in range(n+1):
- for j in range(1,n+1):
- if i**2==1+d*j**2:
- result.append((i,j))
- return result
Utilisation :
- >>> pellfermat(5,2000)
- [(9, 4), (161, 72)]
Et si l’on cherche celui dont la valeur de $y$ est minimale :
- def pellfermat2(d,n):
- result=[]
- for j in range(1,n+1):
- for i in range(1,n+1):
- if i**2==1+d*j**2:
- result.append((i,j))
- return result[0]
Utilisation :
- >>> pellfermat2(7,2000)
- (8, 3)
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 :
- suite = []
- def lehmer1(n):
- if n==0:
- return 12345
- else :
- return (lehmer1(n-1)*40014)%(2**31-85)
- def lehmer2(n):
- if n==0:
- return 67890
- else:
- return (lehmer2(n-1)*40692)%(2**31-249)
- def alea(n):
- alea= (lehmer1(n)-lehmer2(n))%(2**31-86)
- x=alea/(2**31-86)
- suite.append(x)
- return x
Utilisation :
- >>> alea(1)
- 0.9435974024931791
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 :
- from tkinter import *
- def graphe(n):
- #n est le nombre de nombres aléatoires à générer
- fen = Tk()
- taille=200
- can = Canvas(fen, width =taille, height =taille, bg ='ivory')
- can.pack(side =TOP, padx =5, pady =5)
- for i in range(1,n):
- round(alea(i),10)
- #print(round(alea(i),10))
- zoom_x=taille/len(suite)
- zoom_y=taille/max(suite)
- for i in range(n-1):
- 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")
- 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")
Utilisation (Ne cherchez pas à demander trop de points ...) :
- >>> 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 :
- >>> int("AB", 16)
- 171
- >>> int("0xABB", 16)
- 2747
- >>> int("571", 8)
- 377
- >>> int("0571", 8)
- 377
- >>> int("1001010111010",2)
- 4794
- >>> int("0b100101011",2)
- 299
La méthode bin()
convertit et renvoie l’équivalent binaire d’un nombre entier :
- >>> bin(23)
- '0b10111'
La fonction hex()
convertit un nombre entier number en nombre hexadecimal qui lui correspond.
- >>> hex(255)
- '0xff'
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$ :
- def decimalversbase(n,b):
- res=""
- memoire=n
- while memoire>0:
- res=str(memoire%b)+res
- memoire=memoire//b
- return res
- >>> decimalversbase(1000,3)
- '1101001'
- >>> decimalversbase(1000,8)
- '1750'
- >>> decimalversbase(1000,2)
- '1111101000'
Et une autre pour convertir un nombre dans une base $p <10$ en un entier décimal :
- def baseversdecimal(n,b):
- res=0
- if n<b:
- res=n
- else:
- n=str(n)
- l=len(str(n))
- for i in range(l):
- res+=int(n[l-i-1])*b**i
- return res
Utilisation :
- >>> baseversdecimal(1101001,3)
- 1000
- >>> baseversdecimal(1750,8)
- 1000
Et en combinant les deux on obtient la conversion d’un nombre en base $a < 10$ en un nombre en base $b < 10$ :
- def baseaversbaseb(n,a,b):
- return decimalversbase(baseversdecimal(n,a),b)
Utilisation :
- >>> baseaversbaseb(1750,8,3)
- '1101001'
Pour des conversions dans des bases inférieures à 36 :
- from math import pow, log
- def convert(n,a,b):
- n=str(n)
- " n est un str, a et b des entiers "
- 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']
- nb_decimal = 0
- i = len(n)-1
- for x in n:
- nb_decimal += chiffres.index(x)*int(pow(a,i))
- i= i-1
- nbr_chars = int(log(nb_decimal)/log(b)+1)
- res = ""
- for i in range(0, nbr_chars):
- x = int(nb_decimal / pow(b,nbr_chars-1-i))
- nb_decimal -= int(x * pow(b,nbr_chars-1-i))
- res = res+chiffres[x]
- return res
Utilisation :
- >>> convert(1000, 12, 5)
- '23403'
- >>> convert(1000, 35, 5)
- '2333000'
Fermat : « Tout nombre est carré ou somme de 2, 3 ou 4 carrés. »
On peut alors chercher un algorithme qui donne cette décomposition avec le moins de termes possibles :
- from math import sqrt
- def fermat(n):
- N=round(sqrt(n))+1
- for i in range(1,N):
- if i**2==n:
- return str(i)+"^2"
- for i in range(1,N):
- for j in range(1,N):
- if i**2+j**2==n:
- return str(i)+"^2 + "+str(j)+"^2"
- for i in range(1,N):
- for j in range(1,N):
- for k in range(1,N):
- if i**2+j**2+k**2==n:
- return str(i)+"^2 + "+str(j)+"^2 + "+str(k)+"^2"
- else:
- for l in range(1,N):
- if i**2+j**2+k**2+l**2==n:
- return str(i)+"^2 + "+str(j)+"^2 + "+str(k)+"^2 + "+str(l)+"^2"
Utilisation :
- >>> fermat(36)
- '6^2'
- >>> fermat(17)
- '1^2 + 4^2'
- >>> fermat(21)
- '1^2 + 2^2 + 4^2'
- >>> fermat(23)
- '1^2 + 2^2 + 3^2 + 3^2'
- >>> fermat(412)
- '1^2 + 1^2 + 7^2 + 19^2'
- >>> fermat(4123)
- '1^2 + 1^2 + 5^2 + 64^2'
Graphes et matrices
MATRICES
Calculs de base
Nous sommes tentés de faire cela :
- >>> A=[[2,-1,3],[4,-2,0],[-3,3,6]]
- >>> B=[[ 1,3,11],[7,9,17],[13,15,5]]
- >>> 3*A-2*B
- Traceback (most recent call last):
- File "<pyshell#37>", line 1, in <module>
- 3*A-2*B
- TypeError: unsupported operand type(s) for -: 'list' and 'list'
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 :
- >>> A=[0]
- >>> 10*A
- [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Forcément le calcul de 3*A-2*B
allait poser problème ...
Une solution est d’utiliser le module numpy
:
- >>> import numpy as np
- >>> A=[[2,-1,3],[4,-2,0],[-3,3,6]]
- >>> A
- [[2, -1, 3], [4, -2, 0], [-3, 3, 6]]
- >>> A=np.array(A)
- >>> A
- array([[ 2, -1, 3],
- [ 4, -2, 0],
- [-3, 3, 6]])
La commande array()
tranforme la liste de listes A en un tableau A.
On peut maintenant calculer :
- >>> B=np.array([[ 1,3,11],[7,9,17],[13,15,5]])
- >>> 3*A-2*B
- array([[ 4, -9, -13],
- [ -2, -24, -34],
- [-35, -21, 8]])
- >>> -A
- array([[-2, 1, -3],
- [-4, 2, 0],
- [ 3, -3, -6]])
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/A
- Warning (from warnings module):
- File "<pyshell#50>", line 1
- RuntimeWarning: divide by zero encountered in true_divide
- array([[ 0.5 , -1. , 0.33333333],
- [ 0.25 , -0.5 , inf],
- [-0.33333333, 0.33333333, 0.16666667]])
- >>> A**2
- array([[ 4, 1, 9],
- [16, 4, 0],
- [ 9, 9, 36]])
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 :
- >>> np.linalg.inv(A)
- array([[-0.66666667, 0.83333333, 0.33333333],
- [-1.33333333, 1.16666667, 0.66666667],
- [ 0.33333333, -0.16666667, 0. ]])
Si l’on veut calculer le produit matriciel A.B
:
- >>> np.dot(A,B)
- array([[ 34, 42, 20],
- [-10, -6, 10],
- [ 96, 108, 48]])
Autre possibilité :
- >>> A.dot(B)
- array([[ 34, 42, 20],
- [-10, -6, 10],
- [ 96, 108, 48]])
Si l’on veut calculer le cube de A :
- >>> np.linalg.matrix_power(A,3)
- array([[-54, 63, 117],
- [-36, 36, 72],
- [-45, 63, 126]])
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.
$
- >>> A = np.array([[2, 1], [-1, 4]])
- >>> b= np.array([15, -7])
- >>> np.linalg.solve(A, b)
- array([7.44444444, 0.11111111])
On peut bien sûr se passer de cette commande puisque si A est inversible $Ax = b \Leftrightarrow x = A^{-1}b$ :
- >>> np.dot(np.linalg.inv(A),b)
- array([7.44444444, 0.11111111])
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}$
- def inversion2(m11,m12,m21,m22):
- d=m11*m22-m12*m21
- if d==0:
- return "La matrice M n'est pas inversible car son déterminant est nul."
- else:
- return [[m22/d,-m12/d],[-m21/d,m11/d]]
Utilisation :
- >>> inversion2(1,0,-2,3)
- [[1.0, 0.0], [0.6666666666666666, 0.3333333333333333]]
- >>> inversion2(1,-1,-2,3)
- [[3.0, 1.0], [2.0, 1.0]]
- >>> inversion2(1,-1,-3,3)
- "La matrice M n'est pas inversible car son déterminant est nul."
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 :
- >>> np.zeros(3)
- array([0., 0., 0.])
- >>> np.zeros((3,3))
- array([[0., 0., 0.],
- [0., 0., 0.],
- [0., 0., 0.]])
- >>> np.zeros((3,2))
- array([[0., 0.],
- [0., 0.],
- [0., 0.]])
- >>> np.eye(3)
- array([[1., 0., 0.],
- [0., 1., 0.],
- [0., 0., 1.]])
- >>> np.eye(3,2)
- array([[1., 0.],
- [0., 1.],
- [0., 0.]])
- >>> np.ones(3)
- array([1., 1., 1.])
- >>> np.ones((3,2))
- array([[1., 1.],
- [1., 1.],
- [1., 1.]])
- >>> np.diag([-2,3,1])
- array([[-2, 0, 0],
- [ 0, 3, 0],
- [ 0, 0, 1]])
- >>> np.diag([-2,3,1],1)
- array([[ 0, -2, 0, 0],
- [ 0, 0, 3, 0],
- [ 0, 0, 0, 1],
- [ 0, 0, 0, 0]])
- >>> np.diag([-2,3,1],-2)
- array([[ 0, 0, 0, 0, 0],
- [ 0, 0, 0, 0, 0],
- [-2, 0, 0, 0, 0],
- [ 0, 3, 0, 0, 0],
- [ 0, 0, 1, 0, 0]])
- >>> np.random.rand(5)
- array([0.5704163 , 0.46493989, 0.85085904, 0.56177765, 0.06711559])
- >>> np.random.rand(5,2)
- array([[0.69873188, 0.66611399],
- [0.76550507, 0.25193332],
- [0.28569524, 0.67081913],
- [0.63928606, 0.97949063],
- [0.30312998, 0.31526063]])
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 :
- import numpy as np
- def distribution_invariante(lambd,mu):
- A=np.array([[-lambd,mu],[1,1]])
- B=np.array([0,1])
- pi=np.dot(np.linalg.inv(A),B)
- return pi
Utilisation :
- >>> distribution_invariante(0.4,0.7)
- array([0.63636364, 0.36363636])
Vérification (On regarde si $\pi P = \pi$) :
- >>> np.dot([0.63636364, 0.36363636],[[0.6,0.4],[0.7,0.3]])
- array([0.63636364, 0.36363636])
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 :
- from random import randint
- def ehrenfest(n,k):
- L=[1]*k #On met les k boules dans l'urne 1
- for i in range(n):
- choix=randint(0,k-1) #On choisit une boule au hasard
- if L[choix]==1:
- L[choix]=2 #Si elle est dans l'urne 1, on la met dans la 2
- else:
- L[choix]=1 #Si elle est dans l'urne 2, on la met dans la 1
- return L
Utilisation :
- >>> ehrenfest(1000,10)
- [1, 1, 1, 1, 2, 1, 2, 2, 1, 2]
- >>> ehrenfest(1000,15)
- [2, 2, 1, 1, 1, 2, 1, 2, 1, 1, 2, 1, 1, 2, 1]
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.
- def graphique(n,k):
- L=[1]*k
- x=[0]
- y=[k]
- z=[k//2]*(n+1)
- nb=k
- for i in range(n):
- choix=randint(0,k-1) #On choisit une boule au hasard
- x.append(i)
- if L[choix]==1:
- L[choix]=2 #Si elle est dans l'urne 1, on la met dans la 2
- nb=nb-1
- else:
- L[choix]=1 #Si elle est dans l'urne 2, on la met dans la 1
- nb=nb+1
- y.append(nb)
- pl.plot(x,y)
- pl.plot(x,z)
- pl.show()
- return L.count(1)/k
Utilisation :
- >>> graphique(100,4)
- 0.5
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.
- >>> graphique(10000,500)
- 0.52
Construction de graphes
Avec le module Graphviz, pour obtenir le graphe orienté de la chaine de Markov à deux états vu précédemment :
- from graphviz import Digraph
- g = Digraph('G', filename='Markov', engine='dot')
- g.edge("0", "0",label="1 - μ")
- g.edge("1", "0",label="λ")
- g.edge("1", "1",label="1 - λ")
- g.edge("0", "1",label="μ")
- g.format="svg"
- g.render(view=True)
Utilisation :
- from graphviz import Digraph
- g = Digraph('G', filename='Distrib', engine='dot')
- g.edge("E", "E",label="0.3")
- g.edge("A", "E",label="0.4")
- g.edge("A", "A",label="0.6")
- g.edge("E", "A",label="0.7")
- g.format="svg"
- g.render(view=True)
Utilisation :
Pour un graphe non orienté, pondéré :
- from graphviz import Graph
- g = Graph('G', filename='Itineraire', engine='circo')
- g.edge("Montpellier", "Béziers",label="72 km, 57 min")
- g.edge("Montpellier", "Nîmes",label="56 km, 50 min")
- g.edge("Alès", "Nîmes",label="55 km, 38 min")
- g.edge("Ganges", "Alès",label="51 km, 55 min")
- g.edge("Ganges", "Nîmes",label="67 km, 68 min")
- g.edge("Ganges", "Béziers",label="98 km, 89 min")
- g.edge("Ganges", "Montpellier",label="47 km, 54 min")
- g.edge("Alès", "Montpellier",label="80 km, 96 min")
- g.format="svg"
- g.render(view=True)
Utilisation :
On peut obtenir un rendu différent avec le paramètre engine
(Voir paragraphe Roadmap en suivant ce lien) :
- g = Graph('G', filename='Itineraire', engine='dot')
- g = Graph('G', filename='Itineraire', engine='sfdp')
Annexe : Bernard Parisse suggère des améliorations à l’article qui précède.
Bernard Parisse écrit :
"Je suis en train de regarder les algorithmes d’arithmétique de maths expertes présentés par MathémaTICE, ça me parait bien à part deux suggestions que vous pouvez peut-être faire suivre à son auteur.
* Sur l’algorithme amélioré de factorisation d’entiers : je trouve que ce n’est pas une bonne idée de calculer sqrt(n) pour faire le test d’arrêt, parce qu’on passe par du calcul approché. Il vaut mieux utiliser le test i*i<=n qui, lui, est exact (en plus ça évite d’importer math ).
Évidemment dans le domaine d’application de la factorisation naïve, on ne verra pas le problème, sauf si on est très patient !
Voici ce que ça peut donner, sans les optimisations de chercher de 2 en 2 (et sans utiliser divmod) : <https://www-fourier.univ-grenoble-a...>
* Sur le crible amélioré : il est a mon avis trop difficile à lire et du coup on risque de rater l’essentiel, à savoir qu’on gagne en complexité en s’arrêtant lorsque i*i>n. Le reste n’apporte rien en terme de complexité, c’est de la cuisine typique des langages interprétés pour gagner un facteur constant en temps, en réécrivant une boucle interprétée avec une instruction qui fera la boucle en natif. Si on veut vraiment aller plus vite, il faut de toutes façons utiliser un langage compilé, et les méthodes d’optimisation n’ont rien à voir, en fait les méthodes qu’on utilise pour optimiser en langage interprété sont parfois contre-productives en langage compilé (cf. par exemple https://www-fourier.univ-grenoble-alpes.fr/ parisse/irem/algolycee.html#sec156).
En complément, sur la question d’optimisation en langage interprété vs en langage compilé, dans le cas du crible les lignes en cause sont au debut d’eratosthene2
tableau[2::2] = [False]*((n-2)//2 + n%2)
et dans la boucle la ligne (relativement difficile à lire pour quelqu’un qui n’est pas expérimenté en Python) :
tableau[i::i] = [False]*((n-i)//i + int((n-i)%i>0))
dont le rôle est d’éliminer une boucle en langage interprété pour la faire exécuter en natif par l’interpréteur. Mais cette ligne crée une liste qui peut contenir jusqu’à n/2 éléments, on va donc utiliser 50% de mémoire en plus. Quand on optimise en langage compilé, on est très attentif à utiliser le minimum de mémoire (lire et écrire en mémoire est coûteux, encore plus si on quitte le cache du microprocesseur), on évite absolument de créer une liste, c’est donc tout simplement la boucle présentée dans eratosthene(n) qui est utilisée
for j in range(i*2,n,i) :
tableau[j] = 0 "