Mathématiques avec Python et Ruby/Version imprimable1
Une version à jour et éditable de ce livre est disponible sur Wikilivres,
une bibliothèque de livres pédagogiques, à l'URL :
https://fr.wikibooks.org/wiki/Math%C3%A9matiques_avec_Python_et_Ruby
Nombres en Python
Le langage Python est faiblement typé, ce qui veut dire que c'est au moment de la première affectation d'une variable (son instanciation) que Python devine le type de cette variable. Python possède 5 catégories de nombres:
- int pour les nombres entiers;
- long pour les grands entiers;
- fraction pour les fractions (quotients d'entiers par des entiers);
- decimal et float pour les réels (en réalité des nombres décimaux);
- complex pour les nombres complexes.
Nombres entiers en Python
Les nombres entiers ne sont pas les seuls nombres, comme on le verra dans les chapitres suivants. Alors comment fait Python pour savoir qu'un nombre est entier ? Comme le langage est faiblement typé, il doit le deviner. Le critère est simple : Pour qu'un nombre soit entier, il ne doit pas avoir de virgule (représentée dans Python par un point décimal).
Les nombres entiers dans Python
modifierAinsi, si on entre
a=3
print(type(a))
b=3.14
print(type(b))
c=int(b)
print(c)
on constate que Python sait que a est entier, que b ne l'est pas, et que c peut être entier bien qu'obtenu à partir de b (qui est réel).
Certains calculs devant donner un résultat entier ne le font pas toujours en Python. Par exemple, alors que , Python considère ce nombre comme un réel (non entier)!
from math import *
a=sqrt(100)
print(a.is_integer())
Opérations
modifierAddition, soustraction et multiplication
modifierLes trois premières opérations se notent avec les symboles +, - et * comme dans la plupart des langages de programmation. La somme, la différence et le produit de deux entiers (ou plus) sont des entiers (relatifs):
a=5
b=-8
print(a+b)
print(a-b)
print(a*b)
La tortue
modifierOn peut représenter l'addition par des mouvements successifs de la tortue. Pour cela il faut bien entendu importer le module turtle. Après ça pour additionner 55 et 34, on peut faire avancer la tortue successivement de 55 pixels puis 34 pixels, et regarder où elle est:
from turtle import *
forward(55)
forward(34)
print(position())
Pour soustraire le deuxième nombre au lieu de l'additionner, il suffit de faire reculer la tortue au lieu de la faire avancer; et ça marche même avec des nombres négatifs:
from turtle import *
forward(55)
backward(34)
print(position())
Divisions
modifierIl y a deux sortes de divisions d'entiers en Python: Le quotient euclidien, qui est un entier, et le quotient exact, qui est une fraction (pour Python, un réel):
Quotients
modifierDans les anciennes versions de Python, le script suivant
a=3
b=2
print(a/b)
affichait 1 au lieu de 1.5, parce que pour Python, comme a et b sont entiers, leur quotient était logiquement euclidien. Ceci a changé, maintenant le script ci-dessus produit bien 1.5; mais du coup, si par hasard on voulait quand même calculer une division euclidienne avec Python? Et bien dans ce cas il faudrait dédoubler le slash qui code la division:
a=3
b=2
print(a//b)
Ce peut être utile parce que même en divisant 4 par 2, le résultat est un réel non entier, alors que la division euclidienne produit bien le nombre entier 2.
Reste euclidien
modifierLe reste de la division euclidienne de 13 par 8 est 5. Pour le calculer, on utilise l'opérateur infixé %:
a=13
b=8
print(a%b)
Cette opération permet de travailler sur les congruences. Remarque: Si b=0, on a le même message d'erreur que lorsqu'on divise par 0.
Divisibilité
modifierDeux entiers quelconques ont un pgcd. En Python, on l'obtient avec gcd. Mais cette fonction ne se trouve pas par défaut dans Python; elle se trouve dans le fichier fractions.py, dont on peut importer la totalité des objets par from fractions import *:
a=13572468
b=12345678
print(gcd(a,b))
Comme Python est libre, on peut consulter le source de fractions.py, où on trouve ceci:
def gcd(a, b):
while b:
a, b = b, a%b
return a
On reconnaît (d'autant plus aisément que le langage est concis) l'algorithme d'Euclide.
Puissances
modifierBeaucoup de langages de programmation utilisent le chapeau pour représenter les puissances. Pas Python pour qui le chapeau est déjà pris par une autre opération (le ou exclusif bit à bit). Alors c'est l'astérisque de la multiplication qui est utilisé pour les puissances, mais en le dédoublant:
a=4
b=2
print(a**b)
print(b**a)
Remarque: Si l'exposant est négatif ou non entier, le résultat est un réel. En particulier, les deux opérations ci-dessous ont le même effet:
print(100**0.5)
from math import *
print(sqrt(100))
Priorités opératoires
modifierEn Python comme en algèbre, on effectue dans l'ordre
- Les parenthèses
- Les fonctions (comme l'élévation à une puissance)
- Les multiplications et divisions
- Les additions et soustractions.
Ainsi
print(2+3*5)
affiche 17 et non 25: Les opérations ne sont pas effectuées de gauche à droite, mais en suivant les priorités opératoires.
Itérateurs
modifierUn itérateur de Python est une liste d'entiers. Par exemple, la liste 10, 13, 16, 19, 22, 25 (allant de 3 en 3, de 10 jusqu'à 25) s'appelle range(10,26,3). On peut abréger les itérateurs en omettant le troisième argument (par défaut 1) voire le premier (par défaut 0). Ainsi range(5) contient les entiers 0, 1, 2, 3 et 4 (0 compris, 5 non compris, ce sont donc bien les 5 premiers entiers naturels).
Lorsqu'on lance un dé, le résultat obtenu peut être 1, 2, 3, 4, 5 ou 6. Il peut donc être décrit par range(1,7). Pour vérifier que l'évènement "le dé tombe sur 3" est possible alors que l'évènement "le dé tombe sur 7" est impossible, on peut faire
dice=range(1,7)
print(3 in dice)
print(7 in dice)
Fractions en Python
L'écriture de nombres non entiers sous forme de fractions a de loin précédé celle des nombres décimaux, puisque les égyptiens et les babyloniens les utilisaient déjà. Chaque fois que le dénominateur n'est pas une puissance de 10, on continue encore de nos jours à utiliser des écritures fractionnaires plus ou moins cachées, comme dans les exemples suivants :
- Lorsqu'on dit qu'un homme mesure 5 pieds 7 pouces, ça signifie que sa taille, en pieds, est (il y a douze pouces dans un pied);
- Lorsqu'on dit qu'il est 8 heures 13, c'est que depuis minuit, il est passé exactement heures (soit 493 minutes).
- Lorsque Roméo se plaint d'avoir attendu Juliette pendant plus de trois quarts-d'heure, il exprime la durée de son attente insoutenable sous la forme d'une fraction...
- Les probabilités se donnent aussi souvent sous forme de fractions (le plus souvent des fractions égyptiennes). Comme dans "j'ai une chance sur 10 millions de recevoir une météorite sur la tête" ou "Casaque Jaune est donné favori à 5 contre 1".
- Idem parfois pour les statistiques : "5 Français sur 7 estiment qu'il y a trop de sondages sur les sondages"...
L'égalité 0,2+0,5=0,7 peut s'écrire mais l'égalité ne peut pas s'écrire sous forme décimale exacte parce que le résultat n'est pas décimal. Python affiche 1/2+1/3=0.8333333333333333 et malgré l'abondance de chiffres, cette égalité n'est pas exacte.
Pour faire des calculs exacts avec des fractions, Python a un module fractions.py. Pour transformer Python en un langage de programmation spécialisé dans les fractions, il suffit de précéder les scripts de ce chapitre de la ligne
from fractions import *
qui importe le module fractions en entier (pourquoi faire dans le détail?). Une alternative est de chercher le fichier __init__.py, initialement vide, et d'y mettre la ligne précédente.
Obtention d'une fraction
modifierPour entrer la fraction dans Python, on entre Fraction(n,d) après avoir importé le module fractions :
from fractions import *
a=Fraction(24,10)
print(a)
On constate que la fraction a été automatiquement simplifiée par Python au moment de son instanciation.
Si on entre 0 comme dénominateur, la fraction ne se crée pas et on a un message d'erreur : Comme une fraction est un quotient, on ne peut pas diviser par 0.
Une fois qu'une fraction est calculée, on peut obtenir son numérateur et son dénominateur par
from fractions import *
a=Fraction(24,10)
print(a.numerator)
print(a.denominator)
Bien entendu, le numérateur de n'est pas 24...
Pour obtenir la valeur de la fraction (le quotient de son numérateur par son dénominateur), on peut lui additionner 0.0 : La somme d'une fraction et d'un réel, même nul, est un réel.
from fractions import *
a=Fraction(24,10)
print(a+0.0)
Réciproquement, on peut convertir un nombre réel en fraction, mais le résultat est parfois surprenant, si on prend un nombre dont le développement décimal s'arrête alors que le développement binaire ne le fait pas :
from fractions import *
a=Fraction.from_float(1.2)
print(a)
On s'attendait à , on a l'approximation un peu surprenante ...
Opérations
modifierLes opérations sur les fractions se notent comme celles sur les autres nombres, mais en général le résultat est une fraction.
Opérations unaires
modifierOpposé
modifierL'opposée d'une fraction s'obtient en la faisant précéder du signe -. Par exemple, la fraction suivante est positive :
from fractions import *
a=Fraction(2,-3)
print(-a)
Inverse
modifierPour obtenir l'inverse d'une fraction, on divise 1 par celle-ci :
from fractions import *
a=Fraction(5,4)
print(1/a)
Addition
modifierLa somme de deux fractions est une fraction :
from fractions import *
a=Fraction(34,21)
b=Fraction(21,13)
print(a+b)
Soustraction
modifierLa différence de deux fractions est une fraction :
from fractions import *
a=Fraction(34,21)
b=Fraction(21,13)
print(a-b)
Multiplication
modifierLe produit de deux fractions est une fraction :
from fractions import *
a=Fraction(34,21)
b=Fraction(21,13)
print(a*b)
Division
modifierLe quotient de deux fractions est une fraction (à condition que la deuxième ne soit pas nulle) :
from fractions import *
a=Fraction(34,21)
b=Fraction(21,13)
print(a/b)
Le reste euclidien continue à être défini pour des fractions, et le résultat est une fraction :
from fractions import *
a=Fraction(32,7)
b=Fraction(7,2)
print(a%b)
Puissance
modifierSi l'exposant est entier, la puissance d'une fraction est une fraction :
from fractions import *
a=Fraction(3,2)
print(a**12)
print(a**(-1))
Mais si l'exposant est un réel, la puissance n'est pas une fraction mais un réel :
from fractions import *
a=Fraction(9,4)
print(a**0.5)
Algorithmes
modifierRéduite de Farey
modifierLa création de réduite de Farey est aisée avec le module fractions de Python:
from fractions import *
def Farey(a,b):
n=a.numerator+b.numerator
d=a.denominator+b.denominator
return Fraction(n,d)
a=Fraction(5,4)
b=Fraction(1,13)
print(Farey(a,b))
C'est si facile qu'on en vient tout de suite à la question suivante : à quoi ça peut bien servir ?
Si on répond que ça sert à fabriquer un arbre de Stern-Brocot, ça n'éclaire peut-être pas beaucoup, mais disons que ça sert à quelque chose de mathématique...
Fractions égyptiennes
modifierUne fraction égyptienne est définie comme une somme d'inverses d'entiers. En effet les égyptiens avaient la réputation de ne pas utiliser de numérateurs. Toute fraction peut s'écrire comme une somme de fractions égyptiennes, et l'algorithme de Fibonacci permet d'en trouver un exemple à partir de la fraction. Dans la version Python ci-dessous, l'algorithme fournit une liste de fractions, toutes de numérateur 1, dont la somme est une fraction donnée f. Mais au cas où f serait supérieure à 1, on commence la liste par un entier :
from fractions import *
from math import *
def egypt(f):
e=int(f)
f-=e
liste=[e]
while(f.numerator>1):
e=Fraction(1,int(ceil(1/f)))
liste.append(e)
f-=e
liste.append(f)
return liste
a=Fraction(21,13)
print(egypt(a))
Quelques explications sur le bricolage ci-dessus : Le dénominateur d'une fraction égyptienne est choisi entier (bien sûr) et plus grand que l'inverse de la fraction f (pour que l'algorithme converge). Une solution serait de prendre la troncature int de l'inverse de f et ajouter 1. Mais si l'inverse de f est entier, on ne doit pas ajouter 1 (sinon la suite est infinie). Alors on utilise la fonction ceil. Donc
- Il a fallu importer le module math qui contient cette fonction ceil;
- Du coup l'objet ceil(1/f) n'est plus un entier mais un réel, et ne peut plus être le dénominateur d'une fraction (message d'erreur de Python). Alors il faut convertir ce réel (qui est déjà entier, mais Python ne le sait pas) en entier, ce qui se fait par int.
Enfin, Python ne possédant pas de boucle do..while, il faut ajouter la dernière fraction égyptienne à la liste, pour que celle-ci soit complète.
En lançant le script ci-dessus, on apprend que .
Nombres réels en Python
Écriture décimale
modifierSi les fractions paraissent si abstraites, c'est sans doute à cause de l'écriture décimale qui est jugée plus concrète que l'écriture fractionnaire.
Nombres décimaux
modifierUn nombre décimal est un nombre dont le développement décimal s'arrête. Un nombre non décimal (comme les premiers réels connus qui ne soient pas des fractions) a donc une infinité de chiffres, et ne peut donc pas être représenté de façon exacte en machine: On est obligé de travailler sur des approximations. Mais même avec des nombres décimaux, on travaille parfois sur des valeurs approchées: En effet les nombres sont représentés en base 2 dans la mémoire de l'ordinateur, et Python ne gère que les nombres binaires, ce que n'est pas 0,1.
Fractions
modifierUne fraction se reconnaît à ce que son développement décimal est périodique à partir d'un certain rang.
print(1/3)
print(1/9)
print(1/11)
print(1/13)
(ces illustrations nécessitent la version 3.2 de Python. En cas de version inférieure, il faut remplacer les numérateurs par des 1.0 pour que la division se fasse dans ).
Pour mieux voir ces chiffres, on peut en afficher plus:
from decimal import *
print(Decimal(1)/Decimal(3))
print(Decimal(1)/Decimal(9))
print(Decimal(1)/Decimal(11))
print(Decimal(1)/Decimal(13))
Irrationnels
modifierPour "construire" un nombre irrationnel (un réel qui ne soit pas une fraction), on peut donc inventer une suite de chiffres qui ne risque pas d'être répétitive, comme avec la constante de Champernowne ou un autre nombre univers.
Nombres algébriques
modifierLe plus vieil irrationnel connu est sans doute la racine carrée de 2:
print(2**0.5)
#autre version:
from math import *
print(sqrt(2))
D'autres connus depuis longtemps sont le nombre d'or, la racine cubique de 2, etc. Les nombres algébriques sont définis comme solutions d'équations polynomiales, donc leurs valeurs décimales approchées sont calculées avec des méthodes comme la dichotomie, la méthode de Newton etc.
Nombres transcendants
modifierDeux nombres transcendants très connus sont e et :
from math import *
print(e)
print(pi)
Voici une manière de calculer une valeur approchée du nombre de Champernowne en Python:
c='0.'
for n in range(1,40):
c+=str(n)
print(float(c))
Exercice pour un premier avril: Calculer une valeur approchée à 3 décimales près de la constante de Chaitin...
Fonctions réelles
modifierOpérations
modifierLes quatre opérations sont notées +, -, * et / en Python, et leur résultat (sauf si on divise par 0) est un réel. On peut aussi calculer le reste euclidien d'un réel par un réel!
from math import *
angle=100%pi
print(angle)
Le signe - désigne aussi l'opposé d'un réel. Sa valeur absolue est notée abs. Sa racine carrée sqrt.
Pour additionner h au nombre x, on peut écrire x+=h au lieu du classique x=x+h.
Puissances, exponentielles et logarithmes
modifierPuissances et exposants
modifierComme pour les entiers, les puissances se notent avec l'astérisque de la multiplication, dédoublé. L'exemple du haut de la page montre comment ceci permet de calculer sans utiliser le module math.
Logarithmes
modifierLe script ci-dessous calcule et affiche l'image de 0,5 par le logarithme népérien, par le logarithme décimal, par les fonctions réciproques du cosinus hyperbolique, du sinus hyperbolique, de la tangente hyperbolique:
print(log(0.5))
print(log10(0.5))
print(acosh(0.5))
print(asinh(0.5))
print(atanh(0.5))
Exponentielles
modifierPour calculer une puissance de 10, on peut utiliser la notation ** s'il s'agit d'un calcul, mais aussi la notation "E" s'il s'agit de l'exposant entier d'un nombre décimal :
a=1.23456E7
print(a)
b=10**2.5
print(b)
Pour calculer l'exponentielle d'un nombre, on peut utiliser exp:
from math import *
a=e**pi
b=exp(pi)
print(a==b)
print(a-b)
Le script suivant calcule les cosinus, sinus et tangente hyperbolique de 2:
from math import *
print(cosh(2))
print(sinh(2))
print(tanh(2))
Fonctions trigonométriques
modifierOn peut convertir des angles de radians en degrés et vice-versa avec degrees(x) et radians(x). Les fonctions trigonométriques directes se notent cos, sin et tan et ces trois fonctions sont en radians. Les fonctions inverses se notent acos, asin et atan et sont aussi en radians.
On peut calculer la fonction de deux variables avec hypot:
from math import *
a=sqrt(3**2+4**2)
b=hypot(3,4)
print(a==b)
Pour connaître l'angle aigu d'un triangle de côtés x et y, on peut, outre le calcul atan(y/x), faire atan2(x,y). Par exemple, si on veut connaître les angles et l'hypoténuse d'un triangle rectangle de côtés 12 cm et 5 cm, on peut utiliser ce script:
from math import *
a=12
b=5
print(degrees(atan2(a,b)))
print(degrees(atan2(b,a)))
print(hypot(a,b))
Nombres complexes en Python
Python est un langage très utilisé dans le domaine scientifique, comme le montre par exemple le choix de SAGE. Et les sciences, en particulier, font grand usage des nombres complexes, essentiellement depuis leur choix par Cauchy. Les physiciens et les électriciens notant j le nombre complexe dont le carré vaut -1, Python suit ce choix.
Instanciation d'un nombre complexe
modifierDans Python, il suffit d'écrire complex(x,y) pour avoir un nombre complexe :
z=complex(4,3)
print(z)
Même si les coordonnées x et y du point sont entières ou des fractions, elles deviennent des réels lorsque Python instancie le complexe. Voir les propriétés du complexe ci-dessous pour le vérifier.
Si on veut quand même que la lettre i désigne le complexe de carré -1, il suffit de le déclarer comme tel :
i=complex(0,1)
print(i**2)
Opérations
modifierLes quatre opérations se notent respectivement +, -, * et /, et donnent toujours un complexe, même si celui-ci est réel (exemple de la soustraction ci-dessous) :
a=complex(2,3)
b=complex(4,3)
print(a+b)
print(a-b)
print(a*b)
print(a/b)
L'élévation à un exposant se note de la même manière que pour les autres nombres, par **. Mais l'exposant peut même être un complexe!
i=complex(0,1)
print(i**i)
On constate que ...
La racine carrée d'un complexe peut aussi s'obtenir par une élévation de celui-ci à la puissance 0,5 mais dans ce cas on n'obtient qu'une seule des deux racines carrées :
c=complex(7,24)
print(c**0.5)
Mais -4-3i a aussi pour carré 7+24i. Comment fait Python pour choisir entre les deux racines carrées?
Même -1 a deux racines carrées dans , et comme on s'en doute, Python ne choisit pas -i mais i... ou plutôt un complexe proche de celui-ci :
print((-1)**0.5)
Propriétés d'un nombre complexe
modifierLes parties réelle et imaginaire d'un complexe sont des propriétés de l'objet :
z=complex(4,3)
print(z.real)
print(z.imag)
Par contre, le conjugué d'un complexe est une méthode de celui-ci :
z=complex(4,3)
print(z.conjugate())
(on remarque la présence des parenthèses après conjugate)
Forme trigonométrique
modifierPour avoir le module d'un nombre complexe, on entre abs :
z=complex(4,3)
print(abs(z))
Bien entendu, le résultat est réel.
Cependant, pour avoir l'argument de a, il faut charger le module (c'est le cas de le dire!) cmath :
from cmath import *
z=complex(4,3)
print(phase(z))
On remarque que Python utilise le mot phase et non le mot argument. cmath permet aussi de calculer d'un coup le module et l'argument d'un nombre complexe avec polar :
from cmath import *
z=complex(4,3)
print(polar(z))
Pour réaliser l'opération inverse (calculer l'exponentielle d'un nombre imaginaire), on utilise rect :
from cmath import *
print(rect(2,pi/3))
Par exemple, si on veut calculer le plus petit angle et l'hypoténuse d'un triangle rectangle de côtés 12 cm et 5 cm, on peut faire ceci :
from cmath import *
a=12
b=5
z=complex(a,b)
print(phase(z))
print(abs(z))
Fonctions
modifierAvec cmath, on peut appliquer certaines fonctions de la variable réelle à des complexes.
Exponentielles
modifierPour vérifier numériquement que , on peut utiliser l'exponentielle d'un nombre complexe (en l’occurrence, imaginaire) :
from cmath import *
t=complex(0,pi/3)
z=exp(t)
print(z.real==0.5)
print(z.real-0.5)
print(z.imag==sqrt(3)/2)
On voit que la partie réelle n'est pas tout-à-fait égale à 0,5 (la différence est minime mais non nulle), c'est encore une conséquence de la représentation binaire des nombres en machine, puisque le développement binaire de 0,5 est infini, contrairement à son développement décimal.
Le script suivant calcule et affiche les fonctions trigonométriques hyperboliques d'un complexe :
from cmath import *
z=complex(4,3)
print(cosh(z))
print(sinh(z))
print(tanh(z))
Logarithmes
modifierOn peut même calculer le logarithme d'un nombre complexe :
Le script suivant calcule et affiche les fonctions trigonométriques hyperboliques d'un complexe :
from cmath import *
z=complex(4,3)
print(log(z))
Le script suivant calcule et affiche les arguments des fonctions trigonométriques hyperboliques d'un complexe :
from cmath import *
z=complex(4,3)
print(acosh(z))
print(asinh(z))
print(atanh(z))
Fonctions trigonométriques
modifierDirectes
modifierLe script suivant calcule et affiche les fonctions trigonométriques d'un complexe :
from cmath import *
z=complex(4,3)
print(cos(z))
print(sin(z))
print(tan(z))
Inverses
modifierLe script suivant calcule et affiche les arcs des fonctions trigonométriques d'un complexe :
from cmath import *
z=complex(4,3)
print(acos(z))
print(asin(z))
print(atan(z))
Quaternions et octonions en Python
Complexes
modifierOn a vu dans le chapitre précédent que pour Python, un nombre complexe z est essentiellement une structure abritant deux réels, accessibles par z.real et z.imag respectivement. La construction de Cayley-Dickson généralise ce point de vue : En prenant deux complexes a et b, on peut les regrouper dans une nouvelle structure qui est considérée comme un nombre : Un quaternion.
Pour toute la suite, il est conseillé d'importer les fonctions du module math de Python :
from math import *
Mais en réalité, seule la méthode hypot sera utilisée, ce qui signifie que le minimum nécessaire était
from math import hypot
Quaternions
modifierDéfinition
modifierToutes les méthodes (qui permettent de manipuler les quaternions) peuvent être regroupées dans une classe nommée Quaternion :
class Quaternion:
def __init__(self,a,b):
self.a=a
self.b=b
La première méthode, l'initialisation, crée donc deux variables a et b (qui seront des complexes, mais Python ne le sait pas encore) et les rendre accessibles par la notation avec un point, les nombres q.a et q.b (les deux complexes qui définissent le quaternion q) étant des propriétés du quaternion.
Affichage
modifierPour y voir quelque chose, une méthode d'affichage est nécessaire. Comme Python en possède déjà une (la conversion en chaîne de caractères, ou string, notée __str__), on va la surcharger :
def __str__(self):
aff='('
aff+=str(self.a.real)+')+('
aff+=str(self.a.imag)+')i+('
aff+=str(self.b.real)+')j+('
aff+=str(self.b.imag)+')k'
return aff
Un quaternion possède deux propriétés qui sont des nombres complexes, mais chacun d'eux est formé de deux nombres réels, donc un quaternion est formé de 4 nombres réels, ce sont ceux qui sont affichés par la méthode ci-dessus, lorsqu'on entre par exemple print(q).
Fonctions
modifierOpposé
modifier def __neg__(self):
return Quaternion(-self.a,-self.b)
En écrivant -q, on aura désormais l'opposé de q.
Module
modifier def __abs__(self):
return hypot(abs(self.a),abs(self.b))
Grâce à cette méthode, abs(q) retourne un réel, la racine carrée de sa norme.
Conjugué
modifier def conjugate(self):
return Quaternion(self.a.conjugate(),-self.b)
abs(q.conjugate()) renvoie le même résultat que abs(q) parce que tout quaternion a la même norme que son conjugué.
Opérations
modifierAddition
modifierPour additionner deux quaternions, on additionne leurs a respectifs, et leurs b respectifs :
def __add__(self,other):
return Quaternion(self.a+other.a,self.b+other.b)
Soustraction
modifier def __sub__(self,other):
return Quaternion(self.a-other.a,self.b-other.b)
Multiplication
modifierLe produit de deux quaternions est plus difficile à définir :
def __mul__(self,other):
c=self.a*other.a-self.b*other.b.conjugate()
d=self.a*other.b+self.b*other.a.conjugate()
return Quaternion(c,d)
Ce produit admet le quaternion Quaternion(1,0) comme élément neutre, et il est associatif comme tout produit qui se respecte. Mais il n'est pas commutatif :
p=Quaternion(-2+1J,2+3J)
q=Quaternion(3-2J,5+1J)
print(p*q)
print(q*p)
Division
modifierMultiplication par un réel
modifierPour définir le plus facilement possible le quotient de deux quaternions, on a intérêt à définir le produit d'un quaternion par un réel (on peut déjà l'effectuer avec la méthode de produit, en assimilant le réel r avec la quaternion Quaternion(r,0)). Mais comme le symbole de multiplication est déjà utilisé pour la multiplication des quaternions, il en faut un autre pour la multiplication d'un quaternion par un réel. Or il se trouve que __rmul__ (multiplication à l'envers) est encore disponible, donc pour peu qu'on multiplie à droite dans la définition, et à gauche en pratique, on peut ajouter cette méthode :
def __rmul__(self,k):
return Quaternion(self.a*k,self.b*k)
Alors, pour tripler un quaternion q, on peut faire, au choix, q*Quaternion(3,0), ou 3*q.
Division
modifierLe quotient de deux quaternions est désormais simple à définir :
def __div__(self,other):
return self*(1./abs(other)**2*other.conjugate())
Le quotient d'un quaternion par son conjugué est de norme 1 :
p=Quaternion(-2+1J,2+3J)
print(p/p.conjugate())
Cet exemple révèle que , ce qui revient à la décomposition suivante de 81 (un carré) comme somme de 4 carrés : .
Puissances
modifierMême si la multiplication des quaternions n'est pas commutative, elle est associative, c'est tout ce qu'il faut pour définir les puissances des quaternions à exposants entiers (et donc, par formules de Taylor, de la trigonométrie et des exponentielles de quaternions) :
def __pow__(self,n):
r=1
for i in range(n):
r=r*self
return r
Par exemple, on peut calculer le carré d'un quaternion :
q=Quaternion(2J/7,(3+6J)/7)
print(q**2)
Plus généralement, l'ensemble des solutions de est une sphère.
L'étude des itérés de où q et c sont des quaternions, mène à des ensembles de Mandelbrot quaternionniques ([1])
Résumé
modifierVoici le contenu complet de la classe Quaternion de Python :
from math import hypot
class Quaternion:
def __init__(self,a,b):
self.a=a
self.b=b
def __str__(self):
aff='('
aff+=str(self.a.real)+')+('
aff+=str(self.a.imag)+')i+('
aff+=str(self.b.real)+')j+('
aff+=str(self.b.imag)+')k'
return aff
def __neg__(self):
return Quaternion(-self.a,-self.b)
def __add__(self,other):
return Quaternion(self.a+other.a,self.b+other.b)
def __sub__(self,other):
return Quaternion(self.a-other.a,self.b-other.b)
def __mul__(self,other):
c=self.a*other.a-self.b*other.b.conjugate()
d=self.a*other.b+self.b*other.a.conjugate()
return Quaternion(c,d)
def __rmul__(self,k):
return Quaternion(self.a*k,self.b*k)
def __abs__(self):
return hypot(abs(self.a),abs(self.b))
def conjugate(self):
return Quaternion(self.a.conjugate(),-self.b)
def __div__(self,other):
return self*(1./abs(other)**2*other.conjugate())
def __pow__(self,n):
r=1
for i in range(n):
r=r*self
return r
Il suffit de placer tout ça dans un fichier appelé quaternions.py et disposer de toutes ces méthodes en important le module nouvellement créé par
from quaternions import *
Ce qui permet alors de déclarer des quaternions, puis d'effectuer des opérations dessus.
Octonions
modifierCe qui est intéressant avec la construction de Cayley-Dickson utilisée ci-dessus pour les quaternions, c'est qu'elle se généralise : En définissant une structure (un objet) comprenant deux quaternions a et b, on définit un octonion.
Définition et affichage
modifierDéfinition
modifierfrom math import hypot
class Octonion:
def __init__(self,a,b):
self.a=a
self.b=b
Affichage
modifierComme est de dimension 8 sur , l'affichage est plus compliqué que celui des quaternions :
def __str__(self):
aff='('
aff+=str(self.a.a.real)+')+('
aff+=str(self.a.a.imag)+')i+('
aff+=str(self.a.b.real)+')j+('
aff+=str(self.a.b.imag)+')k+('
aff+=str(self.b.a.real)+')l+('
aff+=str(self.b.a.imag)+')li+('
aff+=str(self.b.b.real)+')lj+('
aff+=str(self.b.b.imag)+')lk'
return aff
On voit une arborescence apparaître, le a.a.real désignant la partie réelle du a du quaternion a de l'octonion. La notation avec les points prend ici tout son intérêt, permettant une concision pythonienne qui rendrait presque apprivoisés ces redoutables octonions!
Fonctions
modifierLes fonctions sur les octonions se définissent presque comme celles sur les quaternions, Cayley-Dickson oblige :
Opposé
modifier def __neg__(self):
return Octonion(-self.a,-self.b)
Module
modifier def __abs__(self):
return hypot(abs(self.a),abs(self.b))
C'est pour permettre cette concision qu'on a importé la méthode hypot du module math.
Conjugué
modifier def conjugate(self):
return Octonion(self.a.conjugate(),-self.b)
Opérations
modifierAddition
modifier def __add__(self,other):
return Octonion(self.a+other.a,self.b+other.b)
Encore une fois, le fait d'avoir surchargé la méthode __add__ de Python permet de noter simplement m+n la somme des octonions m et n.
Soustraction
modifier def __sub__(self,other):
return Octonion(self.a-other.a,self.b-other.b)
Multiplication
modifier def __mul__(self,other):
c=self.a*other.a-other.b*self.b.conjugate()
d=self.a.conjugate()*other.b+other.a*self.b
return Octonion(c,d)
Non seulement la multiplication des octonions n'est pas commutative, elle n'est plus associative non plus :
m=Octonion(Quaternion(3+4J,2-7J),Quaternion(1+3J,5-3J))
n=Octonion(Quaternion(2+1J,1-3J),Quaternion(2-2J,1+1J))
o=Octonion(Quaternion(3-2J,-5+3J),Quaternion(1-2J,2-1J))
print((m*n)*o)
print(m*(n*o))
Division
modifierGrâce à la notion de conjugué, on peut facilement définir le quotient de deux octonions, mais c'est encore plus facile en multipliant un octonion par un réel :
Produit par un réel
modifier def __rmul__(self,k):
return Octonion(k*self.a,k*self.b)
Quotient de deux octonions
modifier def __div__(self,other):
return self.conjugate()*(1./abs(other)**2*other)
Là encore, le quotient d'un octonion par son conjugué est de norme 1 :
m=Octonion(Quaternion(3+4J,2-7J),Quaternion(1+3J,5-3J))
n=Octonion(Quaternion(2+1J,1-3J),Quaternion(2-2J,1+1J))
print(m/m.conjugate())
print(abs(n/n.conjugate()))
Ces calculs permettent de décomposer un carré en somme de 8 carrés.
Puissances
modifierComme la multiplication des octonions n'est pas associative, les puissances des octonions ne présentent guère d'intérêt, sauf pour la puissance 2, et sur , l'ensemble des solutions de l'équation est une sphère de dimension 6.
Résumé
modifierLa classe Octonion de Python peut se résumer à ceci :
class Octonion:
def __init__(self,a,b):
self.a=a
self.b=b
def __str__(self):
aff='('
aff+=str(self.a.a.real)+')+('
aff+=str(self.a.a.imag)+')i+('
aff+=str(self.a.b.real)+')j+('
aff+=str(self.a.b.imag)+')k+('
aff+=str(self.b.a.real)+')l+('
aff+=str(self.b.a.imag)+')li+('
aff+=str(self.b.b.real)+')lj+('
aff+=str(self.b.b.imag)+')lk'
return aff
def __neg__(self):
return Octonion(-self.a,-self.b)
def __add__(self,other):
return Octonion(self.a+other.a,self.b+other.b)
def __sub__(self,other):
return Octonion(self.a-other.a,self.b-other.b)
def __mul__(self,other):
c=self.a*other.a-other.b*self.b.conjugate()
d=self.a.conjugate()*other.b+other.a*self.b
return Octonion(c,d)
def __rmul__(self,k):
return Octonion(k*self.a,k*self.b)
def __abs__(self):
return hypot(abs(self.a),abs(self.b))
def conjugate(self):
return Octonion(self.a.conjugate(),-self.b)
def __div__(self,other):
return self.conjugate()*(1./abs(other)**2*other)
Pour peu qu'on l'ait enregistrée dans un fichier octonions.py, il suffit pour pouvoir effectuer des calculs sur les octonions, d'importer ce fichier par
from octonions import *
Bibliographie
modifier- De par leur utilité en infographie 3D, les quaternions sont utilisés dans Blender, avec cette description : [2]
- Sur les octonions, le livre de John Baez est une lecture hautement conseillée : [3]
Python et probabilités
Ensembles en Python
Dans la théorie des probabilités telle qu'elle a été axiomatisée par Kolmogorov, un évènement est noté par la liste des éventualités qui le réalisent, notée entre accolades.
Représentation des évènements en Python
modifierÉvènements certain et impossible
modifierL'évènement impossible est noté .
L'évènement certain (noté ) est la liste de toutes les éventualités. Pour savoir si un élément est dans un ensemble, on utilise le mot-clé in comme dans 6 in omega qui est un booléen.
Avec un dé
modifierOn s'apprête à lancer un dé. Alors l'évènement "le résultat sera plus petit que 5" est décrit par l'ensemble . De même, l'évènement "le résultat sera pair" est représenté par .
En Python cela donne :
univers={1,2,3,4,5,6}
petit={1,2,3,4}
pair={2,4,6}
Avec des cartes
modifierCette fois-ci, on extrait au hasard une carte parmi un jeu de 32 cartes.
Faire la liste des 32 cartes (pour constituer l'univers) est un peu fastidieux, alors on va laisser Python le faire :
valeurs={1,7,8,9,10,'Valet','Dame','Roi'}
couleurs={'carreau','cœur','pique','trèfle'}
univers={0}
for v in valeurs:
for c in couleurs:
univers.add(str(v)+' '+c)
univers.remove(0)
print(len(univers))
Il a été nécessaire de mettre initialement un 0 dans l'univers, puis de l'enlever à la fin. C'est pour tromper le typage faible de Python qui considère les accolades vides comme un objet de type dictionnaire et non comme un ensemble. De plus, on transforme les valeurs des cartes en texte même si ce sont des nombres.
L'évènement "la carte est une figure" (pas un nombre) se construit par
couleurs={'carreau','cœur','pique','trèfle'}
figure={0}
for v in {'Valet','Dame','Roi'}:
for c in couleurs:
figure.add(v+' '+c)
figure.remove(0)
print(univers)
Et l'évènement "la carte est un pique" se construit de manière analogue :
valeurs={1,7,8,9,10,'Valet','Dame','Roi'}
pique={0}
for v in valeurs:
pique.add(str(v)+' pique')
pique.remove(0)
print(pique)
Calcul d'évènements
modifierÉvènements simultanés
modifierNotation
modifierL'évènement "A et B" se note , et l'opération se note en Python par une esperluette (&) qui est d'ailleurs une ancienne représentation du mot et en latin.
Avec le dé
modifierunivers={1,2,3,4,5,6}
petit={1,2,3,4}
pair={2,4,6}
print(petit&pair)
Avec les cartes
modifierprint(figure&pique)
L'affichage confirme qu'il n'y a que trois cartes qui sont à la fois des figures et des piques : Les trois figures de pique Ogier, Pallas et David.
Disjonction
modifierNotation
modifierDe même l'évènement "A ou B" se note , et en Python, le symbole pipe (trait vertical). Python enlève automatiquement les doublons.
Avec le dé
modifierunivers={1,2,3,4,5,6}
petit={1,2,3,4}
pair={2,4,6}
print(petit|pair)
Avec les cartes
modifierprint(figure|pique)
print(len(figure|pique))
On peut compter les 17 cartes à la main, mais le comptage par Python est plus sûr. On constate que pour Python, le nombre d'éventualités d'un évènement s'appelle sa longueur.
Contraire d'un évènement
modifierPour calculer le contraire d'un évènement, on le soustrait à l'univers.
Avec le dé
modifierunivers={1,2,3,4,5,6}
petit={1,2,3,4}
pair={2,4,6}
print(univers-petit)
print(univers-pair)
On constate que le contraire de pair est impair...
Avec les cartes
modifierprint(univers-figure)
print(univers-pique)
Probabilités
modifierOn a vu ci-dessus que pour Python, le nombre d'éléments d'un évènement est appelé sa longueur. On peut alors définir la probabilité d'un évènement comme le quotient de sa longueur par celle de l'univers.
Avec le dé
modifierdef proba(evenement):
return len(evenement)/len(univers)
print(proba(petit))
print(proba(pair))
print(proba(univers/pair&petit))
p1=proba(petit)+proba(pair)-proba(petit&pair)
p2=proba(petit|pair)
print(p1==p2)
Avec les cartes
modifierdef proba(evenement):
return len(evenement)/len(univers)
print(proba(figure))
print(proba(pique))
print(proba(univers))
p1=proba(figure)+proba(pique)-proba(figure&pique)
p2=proba(figure|pique)
print(p1==p2)
Nombres pseudoaléatoires en Python
Pour faire des simulations en proba, on a besoin de nombres pseudo-aléatoires.
Obtention
modifierRecette
modifierVersions jusqu'à la 2.7
modifierPython 2.7 calculait les nombres pseudo-aléatoires avec l'algorithme de Wichmann-Hill, basé sur 3 générateurs congruentiels linéaires en parallèle (extrait du code source - on remarquera l'allusion à Fermat dans l'avant-dernière ligne):
# This part is thread-unsafe:
# BEGIN CRITICAL SECTION
x, y, z = self._seed
x = (171 * x) % 30269
y = (172 * y) % 30307
z = (170 * z) % 30323
self._seed = x, y, z
# END CRITICAL SECTION
# Note: on a platform using IEEE-754 double arithmetic, this can
# never return 0.0 (asserted by Tim; proof too long for a comment).
return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0
Ceci suggère quelques exercices d'arithmétique:
- Vérifier que les modules 30269, 30307 et 30323 sont premiers.
- Vérifier que les multiplicateurs 171, 172 et 170 sont primitifs modulo leurs modules respectifs.
- Simuler l'un des générateurs congruentiels (x par exemple est une suite géométrique de raison 171 dans le corps de Galois ).
- Il résulte du point 2 ci-dessus que la suite x est périodique de période 30268. L'addition de la dernière ligne donnerait une période de . Calculer ce nombre et le décomposer en facteurs premiers...
Version 3.2
modifierPython 3.2 utilise le Mersenne twister, implémenté en C. La période est annoncée comme égale à . Démontrer la primalité de ce nombre de Mersenne est difficile, mais un exercice intéressant est de démontrer la primalité de son exposant 19 937.
Une fois qu'on sait construire un nombre entier pseudo-aléatoire compris entre 0 et un entier fixe N-1, il suffit de le diviser par N pour simuler un réel pseudo-aléatoire compris entre 0 (inclus en théorie) et 1 (exclu). Ce réel pseudo-aléatoire est celui appelé random et sa loi est uniforme sur l'intervalle allant de 0 à 1. Une transformation affine permet alors d'obtenir à partir de lui un réel pseudo-aléatoire uniforme entre deux réels donnés, puis une troncature permet de retrouver des entiers aléatoires, modulo un nombre entier donné.
Le module random de Python 3.2 ne sait pas simuler de variables aléatoires binomiales ou de Poisson, mais il fournit des variables aléatoires normales avec l'algorithme de Kinderman et Monahan :
NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)
def normalvariate(self, mu, sigma):
random = self.random
while 1:
u1 = random()
u2 = 1.0 - random()
z = NV_MAGICCONST*(u1-0.5)/u2
zz = z*z/4.0
if zz <= -_log(u2):
break
return mu + z*sigma
Il est basé sur deux variables aléatoires u1 et u2 uniformes entre 0 et 1 (sauf que u2 ne peut être nulle). u1 est ensuite centrée, et divisée par u2. Le tout est normalisé avec la constante magique , élevé au carré et on finit par lui appliquer une transformation affine pour que ses espérance et écart-type soient ceux désirés.
En plus, une variable aléatoire normale centrée réduite peut être obtenue avec gauss qui utilise l'algorithme de Box-Muller.
Python 3.2 possède aussi des simulateurs de tirage sans remise, et même de permutations aléatoires, où on mélange les cartes en multipliant les échanges de deux cartes (technique analogue à celle consistant à couper le jeu plusieurs fois de suite):
for i in reversed(range(1, len(x))):
# pick an element in x[:i+1] with which to exchange x[i]
j = int(random() * (i+1))
x[i], x[j] = x[j], x[i]
Cet algorithme illustre le fait que toute permutation est engendrée par des transpositions.
Syntaxe
modifierToutes ces merveilles se trouvant dans le module random, on doit charger celui-ci pour faire des simulations avec Python:
from random import *
Variables uniformes
modifierContinue
modifierLe plus simple à obtenir c'est le fameux nombre pseudo-aléatoire entre 0 et 1:
from random import *
print(random())
Pour avoir un nombre entre -2 et 3, on peut faire
from random import *
print(uniform(-2,3))
Entière
modifierPour lancer un dé, on peut utiliser l'une des méthodes suivantes:
from random import *
print(randint(1,6))
print(randrange(1,7))
print(choice(range(1,6)))
from math import *
print(ceil(6*random()))
Variables continues
modifierLa loi exponentielle de paramètre quelconque (ici 3) peut être simulée:
from random import *
print(expovariate(3))
Une variable gaussienne peut être simulée par deux méthodes (ici centrée et réduite):
from random import *
print(normalvariate(0,1))
print(gauss(0,1))
random simule aussi d'autres lois continues comme les lois Beta, Gamma, de Weibull etc.
Variables binomiales et de Poisson
modifierOn l'a vu plus haut, random ne simule pas de variables binomiales. Mais une variable binomiale de paramètres n et p pouvant être définie comme somme de n variables de Bernoulli de probabilité p indépendantes entre elles, est simulable par boucle:
from random import *
def Bernoulli(p):
if random()<p:
return 1
else:
return 0
def binomial(n,p):
somme=0
for k in range(n):
somme+=Bernoulli(p)
return somme
print(binomial(25,0.2))
Cet algorithme nécessite d'être amélioré, en effet rien n'empêche de l'appeler avec des valeurs de p supérieures à 1 ou des valeurs de n non entières.
Pour simuler une loi de Poisson, on peut utiliser une variable binomiale qui l'approche, comme par exemple
def Poisson(l):
return binomial(1000000,l/1000000)
mais c'est très long à calculer. Une méthode plus rapide utilise la loi exponentielle de même paramètre:
from random import *
def Poisson(l):
t=1.0/expovariate(l)
return int(t)
Permutations et arrangements
modifierPermutations
modifierEn 1741, dans les Mémoires de l'Académie des Sciences de Berlin, Leonhard Euler publiait un article titré Calcul de la probabilité du jeu de rencontre. Le nom initial du jeu de rencontre était jeu de treize parce qu'il se jouait à 13 cartes. Mais Euler généralise le jeu à n cartes.
Le jeu de rencontre est un jeu de hasard, où deux personnes ayant chacune un entier jeu de cartes, en tirent à la fois une carte après l'autre, jusqu'à ce qu'il arrive, qu'elles rencontrent la même carte: et alors l'une des deux personnes gagne. Or, lorsqu'une telle rencontre n'arrive point du tout, alors c'est l'autre des deux personnes qui gagne. Cela posé, on demande la probabilité, que l'une et l'autre de ces deux personnes aura de gagner.
Dans cet excellent article (16 pages en Français),
Pourvu donc que le nombre de cartes ne soit pas moindre que 12, l'espérance de A sera toujours à celle de B à peu près comme 12 à 7 ... Ou bien parmi 19 jeux qu'on joue, il y en aura probablement 12 qui font gagner A, et 7 qui feront gagner B.
On constate qu'Euler utilisait le mot espérance là où aujourd'hui on écrit probabilité, et sa conclusion s'écrit algébriquement et , expliquant pourquoi dans les diligences de l'époque, les gens pariaient à 12 contre 7 sur une rencontre dans ce jeu.
Pour vérifier cela avec 2 jeux de 32 cartes, ça peut prendre du temps (bien que le jeu s'arrête dès la première rencontre). Alors la simulation offerte par Python permet de rapidement simuler un grand nombre de parties (Python est une sorte de diligence supersonique).
Pour mélanger un jeu de 32 cartes, on doit d'abord le construire, avec la technique vue précédemment, sauf que le mélange ne peut être fait que sur une liste et pas sur un ensemble:
valeurs={1,7,8,9,10,'Valet','Dame','Roi'}
couleurs={'carreau','cœur','pique','trèfle'}
jeu1=[str(v)+' '+c for v in valeurs for c in couleurs]
jeu2=[str(v)+' '+c for v in valeurs for c in couleurs]
from random import *
shuffle(jeu2)
diagnostic='pas de rencontre'
for i in range(32):
if jeu1[i]==jeu2[i]:
diagnostic='rencontre sur la '+str(i)+'ème carte, qui est le '+jeu1[i]
break
print(diagnostic)
On a créé deux jeux de 32 cartes, mélangé l'un des deux (jeu2) puis comparé carte à carte les deux jeux. Par défaut le texte de sortie est pas de rencontre. En effet c'est seulement s'il y a une rencontre qu'il est modifié, et remplacé par le nom et le numéro de la carte commune aux deux jeux. L'instruction break sert à éviter de continuer à jouer après la fin du jeu.
Tirage sans remise
modifierUn magicien demande à une personne de l'assistance de tirer une carte d'un jeu de 32. Quelle est la probabilité que ce soit l'as de pique?
Pour simuler l'expérience, on va utiliser la fonction choice qui permet de tirer une carte. On va alors construire le jeu de cartes comme dans l'article précédent, à ceci près que choice attend une liste et non un ensemble:
valeurs={1,7,8,9,10,'Valet','Dame','Roi'}
couleurs={'carreau','cœur','pique','trèfle'}
univers=[str(v)+' '+c for v in valeurs for c in couleurs]
from random import *
for n in range(100):
if choice(univers)=="1 pique":
print('victoire !')
Une estimation de la probabilité de l'évènement peut alors se faire en comptant le nombre de fois que le mot victoire a été écrit: Il est, en pourcents, la fréquence de l'évènement.
Pour jouer au poker, on peut simuler le choix d'une main par un tirage de 5 éléments (sans répétition) parmi les 32:
valeurs={1,7,8,9,10,'Valet','Dame','Roi'}
couleurs={'carreau','cœur','pique','trèfle'}
univers=[str(v)+' '+c for v in valeurs for c in couleurs]
from random import *
hand=sample(univers,5)
print(hand)
Simulation avec Python
Avec des nombres pseudo-aléatoires et des boucles qui permettent de répéter un grand nombre de fois une expérience élémentaire, on peut simuler des phénomènes aléatoires, et utiliser la loi des grands nombres pour estimer des probabilités (parfois difficiles voire impossibles à calculer).
Pour faire des statistiques, on a besoin de tableaux (d'effectifs). Une difficulté supplémentaire apparaît alors: la première valeur d'un tableau est donnée par l'indice 0 et non l'indice 1 de celui-ci. Ce qui oblige parfois à des décalages d'indice.
Lancers de dés
modifierUn dé
modifierPour vérifier que le dé virtuel défini par Python est équilibré, on peut le lancer un grand nombre de fois (par exemple 6000) et compter les différents résultats:
from random import *
effectifs=[0,0,0,0,0,0]
for n in range(6000):
dice=randint(1,6)
effectifs[dice-1]+=1
print(effectifs)
On peut simplifier la création du tableau en multipliant par 6 un tableau d'une seule case:
from random import *
effectifs=[0]*6
for n in range(6000):
dice=randint(1,6)
effectifs[dice-1]+=1
print(effectifs)
Deux dés
modifierPour étudier la somme des résultats donnés par deux dés indépendants l'un de l'autre (voir par exemple si elle est équidistribuée), on fait comme avec un seul dé sauf qu'on a deux dés, et qu'on les additionne:
from random import *
effectifs=[0]*11
for n in range(6000):
de1=randint(1,6)
de2=randint(1,6)
twodice=de1+de2
effectifs[twodice-2]+=1
print(effectifs)
Avec des cartes
modifierUne carte
modifierOn tire une carte d'un jeu de 32. Pour estimer la probabilité que ce soit l'as de pique, on répète 3200 fois l'expérience, et on divise le nombre de parties gagnées par 3200:
valeurs={1,7,8,9,10,'Valet','Dame','Roi'}
couleurs={'carreau','cœur','pique','trèfle'}
univers=[str(v)+' '+c for v in valeurs for c in couleurs]
from random import *
somme=0
for n in range(3200):
if choice(univers)=="1 pique":
somme+=1
print(somme/3200.)
print(1./32)
Une main
modifierOn tire 5 cartes d'un jeu de 32. Quelle est la probabilité des évènements suivants:
- On a une couleur (les 5 cartes sont de la même couleur);
- On a un carré d'as (4 des 5 cartes sont des as)?
On reconnaît la couleur d'une carte en regardant les deux dernières lettres de son nom:
valeurs={1,7,8,9,10,'Valet','Dame','Roi'}
couleurs={'carreau','cœur','pique','trèfle'}
univers=[str(v)+' '+c for v in valeurs for c in couleurs]
from random import *
somme=0
for n in range(1000000):
main=sample(univers,5)
couleurs_dans_main={0}
for carte in main:
couleurs_dans_main.add(carte[-2:])
if len(couleurs_dans_main)==1:
somme+=1
print(somme/1000000.)
Les couleurs sont très rares!
On reconnaît un as à ce que son nom commence par un 1 non suivi par un 0 (sinon ce serait un 10). On compte les as de chaque main, et on compte combien de fois on en a 4 (un carré):
valeurs={1,7,8,9,10,'Valet','Dame','Roi'}
couleurs={'carreau','cœur','pique','trèfle'}
univers=[str(v)+' '+c for v in valeurs for c in couleurs]
from random import *
somme=0
for n in range(10000):
main=sample(univers,5)
NombreAs=len([carte for carte in main if carte[0:2]=='1 '])
if NombreAs==4:
somme+=1
print(somme/10000.)
Jeu de rencontre
modifierOn cherche à estimer expérimentalement la probabilité d'une "rencontre" avec deux jeux de 32 cartes (qu'à un moment donné, les deux joueurs, dont l'un a mélangé son jeu, déposent la même carte sur la table). Pour cela, on répète 19000 fois le jeu de rencontre, et on compte combien de rencontres on a eu:
valeurs={1,7,8,9,10,'Valet','Dame','Roi'}
couleurs={'carreau','cœur','pique','trèfle'}
jeu1=[str(v)+' '+c for v in valeurs for c in couleurs]
jeu2=[str(v)+' '+c for v in valeurs for c in couleurs]
from random import *
rencontres=0
for n in range(19000):
shuffle(jeu2)
for i in range(32):
if jeu1[i]==jeu2[i]:
rencontres+=1
break
print(rencontres/19000.)
print(12./19)
Méthode de Monte-Carlo
modifierPour calculer par la méthode de Monte-Carlo, on "crée" un nuage de points à coordonnées uniformes entre 0 et 1, et on compte combien d'entre eux sont à une distance de l'origine inférieure à l'unité. La fréquence de ces points converge vers :
from math import hypot
from random import random
p=len([n for n in range(1000000) if hypot(random(),random())<1])
print(p/1000000*4)
Heureusement, il y a des moyens plus rapides pour calculer !
Statistique inférentielle avec Python
En statistique inférentielle, on cherche à connaître l'inconnu. Pour ce faire, on émet des hypothèses ou on estime des grandeurs partiellement inconnues,
- On se fixe des probabilités a priori de se tromper dans ses estimations ou son test;
- On cherche (difficile) pour quels choix des paramètres (intervalle ou estimateur) ces probabilités sont atteintes;
- On prend ses décisions à partir des choix faits ci-dessus.
La simulation permet de prendre le problème à l'envers, en créant un modèle qui se comporte comme l'échantillon qu'on observe, et en estimant les probabilités difficiles à calculer par des fréquences. Il arrive que, de ces simulations, on puisse inférer (on est là pour ça) des conjectures, qui serviront à élaborer les algorithmes de test ou d'estimation. Certes, dire on prend 2 parce que la simulation a suggéré que, pour 2, la probabilité est de 0,95, ce n'est pas très mathématique, mais la théorie qui disait elle aussi on prend 1,96; concrètement 2 est très souvent basée sur des approximations normales de lois compliquées. On va donc expérimenter et observer sans expliquer ce qu'on observe, mais en se servant de ces observations pour expliquer pourquoi les statisticiens font comme ci et pas comme ça.
Estimations
modifierPour faire des statistiques, il faut un échantillon de données aléatoires ou non. Et pour avoir des données sous Python, le plus simple est de les fabriquer sous Python. Par exemple si on veut faire des statistiques sur les 100 premiers carrés d'entiers, on peut fabriquer une liste contenant ces 100 nombres:
donnees=[n**2 for n in range(100)]
print(len(donnees))
Moyenne
modifierPour calculer la moyenne des nombres qui sont dans donnees, on les additionne et on divise la somme par le nombre de nombres qu'il y a dans donnees:
def moyenne(tableau):
return sum(tableau, 0.0) / len(tableau)
print(moyenne(donnees))
L'algorithme est améliorable puisque si une donnée n'est pas numérique, il ne donne qu'un message d'erreur.
Variance
modifierLa variance est définie comme la moyenne des carrés des écarts à la moyenne:
def variance(tableau):
m=moyenne(tableau)
return moyenne([(x-m)**2 for x in tableau])
print(variance(donnees))
Une variante pour la variance est donnée par la formule de Huyghens: Moyenne des carrés moins le carré de la moyenne.
Écart-type
modifierL'écart-type est défini comme la racine carrée de la variance:
def ecartype(tableau):
return variance(tableau)**0.5
print(ecartype(donnees))
Échantillons
modifierOn peut créer un échantillon de 100 nombres gaussiens d'espérance 16 et d'écart-type 2, puis calculer sa moyenne et son écart-type:
from random import *
echantillon=[gauss(16,2) for n in range(100)]
print(moyenne(echantillon))
print(ecartype(echantillon))
On voit que la moyenne est proche de 16 et l'écart-type proche de 2. C'est rassurant. Mais si on y regarde de plus près, on voit un problème: En prenant des échantillons plus petits, on s'attend à ce que leurs moyenne et écart-type fluctuent mais que la moyenne des moyennes (sur beaucoup de petits échantillons) soit 16 et que la moyenne des écarts-types soit proche de 2. C'est vrai pour la moyenne des moyennes mais visiblement pas pour la moyenne des écarts-types:
from random import *
m=[] #liste des moyennes des echantillons
s=[] #liste des ecarts-types des echantillons
for n in range(10000):
echantillon=[gauss(16,2) for k in range(5)]
m.append(moyenne(echantillon))
s.append(ecartype(echantillon))
print(moyenne(m)) # Voisin de 16, c'est rassurant!
print(moyenne(s)) # Largement plus petit que 2!
print(moyenne(s)*2**0.25)
print(ecartype(m)) # Le moyennage resserre les ecarts-types
print(2/5**0.5) # en les divisant par la racine de la taille de l'echantillon
En théorie, le nombre par lequel on doit multiplier la moyenne des écarts-types pour estimer l'écart-type de la population est . Ce n'est pas le cas ici: Il semble que l'algorithme de Python pour effectuer des tirages avec remise introduise un biais. Par contre on découvre expérimentalement ici que la moyenne des écarts-types doit être multipliée par pour avoir un estimateur sans biais de l'écart-type...
Intervalles de confiance
modifierPour des fréquences
modifier
Comme Python sait faire des tirages sans remise, on peut constituer une liste ordonnée de pour et de contre, et y puiser des échantillons au hasard. On peut estimer la proportion d'entre eux qui donne de faux espoirs au maire (au moins 50 pour parmi les 100).
from random import *
population=['contre' for n in range(57000)]+['pour' for n in range(43000)]
shuffle(population)
print(len(population))
#On choisit 1000 echantillons de 100 et on compte combien sont favorables au maire:
p=len([n for n in range(1000) if len([v for v in sample(population,100) if v=='pour'])>=50])
print(p)
print(p/1000)
L'échantillon choisi par l'institut de sondage, s'il a réellement été choisi au hasard, était favorable au maire avec une probabilité d'environ 0,1. Cette probabilité n'est pas si ridicule que ça, et l'institut de sondage aurait pu répondre au maire "c'est la faute à pas de chance": Il est tombé sur les 10 % d'échantillons favorables... Un épisode analogue s'est déroulé lors des élections présidentielles de 1995, le ministre du budget de l'époque ayant un peu trop vite cru aux sondages !
Plus sérieux (et plus prudent, pour éviter la vindicte de l'ancien maire, désormais dans l'opposition, et qui a maintenant le temps de mener une croisade contre les instituts de sondage) eût été la publication par l'institut de sondage, d'un intervalle de confiance, par exemple à 95% (c'est-à-dire un intervalle qui contient en moyenne 95% des échantillons). Expérimentalement, on peut s'inventer un intervalle et compter la fréquence des échantillons de 100 personnes qui sont dedans. Ce sera un estimateur de la probabilité que l'échantillon soit représentatif de l'ensemble de la population:
from random import *
h=0.1
p=0
for n in range(1000):
pourcentage=len([v for v in sample(population,100) if v=='pour'])/100
if pourcentage>0.43-h and pourcentage<0.43+h:
p+=1
print(p/1000)
On voit que l'intervalle [0,33 ; 0,53] obtenu avec h=0,1 est un intervalle à 95 %. En modifiant la valeur de h on constate que si h diminue (l'intervalle rétrécit), on perd de la confiance (la probabilité qu'il soit bon diminue aussi). On trouve par tâtonnements la valeur de h pour laquelle la confiance de l'intervalle vaut 95 %, puis par changement de la taille de l'échantillon, on peut conjecturer le lien entre h et la taille de l'échantillon.
Pour des moyennes
modifierLà encore, on peut facilement tester des intervalles de confiance, pour des moyennes de variables aléatoires normales, par exemple d'espérance 16, d'écart-type 2 et indépendantes entre elles:
from random import *
h=0.1
p=0
for n in range(1000):
echantillon=[]
for k in range(100):
echantillon.append(gauss(16,2))
m=moyenne(echantillon)
if m>16-h and m<16+h:
p+=1
print(p/1000)
On découvre que l'intervalle de confiance [15,9 ; 16,1] donné ci-dessus (pour h=0,1) est à environ 40% de confiance. En modifiant la valeur de h, on retrouve expérimentalement que pour celle-ci égale à environ , l'intervalle est à 95 % de confiance.
Test d'équirépartition
modifierEn lançant un dé 100 fois, on constate que le 6 est sorti un peu souvent par rapport aux autres nombres:
Résultat du tirage | 1 | 2 | 3 | 4 | 5 | 6 |
Effectifs | 15 | 16 | 17 | 16 | 16 | 20 |
On se demande si le dé est équilibré. Pour cela, on se choisit comme critère de test la somme des carrés des écarts aux fréquences théoriques: :
d2=(0.15-1/6)**2+(0.16-1/6)**2+(0.17-1/6)**2+(0.16-1/6)**2+(0.16-1/6)**2+(0.2-1/6)**2
Soit mais encore, que faire avec d2: Est-ce qu'on doit dire que 0,0015 est anormalement élevé et que le dé est truqué, ou que 0,0015 est suffisamment petit pour attribuer ces résultats à la fluctuation d'échantillonnage? Pour le savoir, on va simuler 10000 lancers de dés et calculer l'équivalent de d2 pour chacun d'entre eux, puis faire une étude statistique sur le d2 observé. La réponse statistique à la question Qu'est-ce qui est normal? est en général fournie par les déciles: On dira que le dé est vraisemblablement truqué si le d2 observé est supérieur au neuvième décile de la série. Pour calculer ce décile, on va devoir trier les données.
from random import *
khi2=[]
for n in range(10000):
effectifs=[0]*6
for k in range(100):
effectifs[randint(0,5)]+=1
dsquare=0
for e in effectifs:
dsquare+=(e/100-1/6)**2
khi2.append(dsquare)
khi2.sort()
decile9=khi2[int(0.9*len(khi2))]
print(decile9)
print(d2)
d2 est environ 10 fois plus petit que le neuvième décile de la série, donc on se trompe de plus de 10 % en considérant que le dé est truqué: Il est parfaitement normal pour autant qu'on sache.
Problème des rencontres d'Euler
modifierEuler affirme que, si on joue 19 fois au jeu de treize, le joueur A gagnera probablement 12 fois. On peut détailler cette affirmation, en jouant 19 fois au jeu pour chaque échantillon, et en répétant l'expérience 1000 fois (donc 1000 échantillons de 19 parties). Voici le script:
valeurs={1,7,8,9,10,'Valet','Dame','Roi'}
couleurs={'carreau','cœur','pique','trèfle'}
jeu1=[str(v)+' '+c for v in valeurs for c in couleurs]
jeu2=[str(v)+' '+c for v in valeurs for c in couleurs]
from random import *
pg=[0]*20
for n in range(1000):
rencontres=0
for p in range(19):
shuffle(jeu2)
for i in range(32):
if jeu1[i]==jeu2[i]:
rencontres+=1
break
pg[rencontres]+=1
print(pg)
Le tableau affiché par ce script est celui des parties gagnées par A parmi les 19; il s'agit d'un tableau d'effectifs (le total des nombres entiers affichés est d'ailleurs 1000). Voici un exemple de diagramme en bâtons de ce tableau, traîtreusement déguisé en histogramme:
On voit que le mode de cette série est 12, c'est peut-être ce que voulait dire Euler. Mais peut-être aussi voulait-il dire que l'espérance de la série (en fait, sa moyenne puisqu'on fait de la statistique) est proche de 12. C'est bien le cas, cette moyenne étant proche de 11,9 (avec un écart-type proche de 2). En fait le phénomène aléatoire simulé ci-dessus est assez bien modélisé par une variable aléatoire binomiale de paramètres 19 et , dont l'espérance est 12 et l'écart-type .
Python et analyse
Suites en Python
Une suite est une fonction de dans . Tout ce qui est dit dans le chapitre suivant peut donc être appliqué aux suites. On va donc essentiellement parler de suites récurrentes ici. Numériquement, les suites servent surtout à faire des calculs, selon la méthode suivante:
- Pour calculer un nombre , on invente une suite telle que ;
- On estime que par exemple 1000 est suffisamment proche de pour que .
- On calcule alors (en général par récurrence).
- On considère donc qu'on a fini de calculer .
C'est ainsi par exemple, qu'on calcule
- des racines carrées par la méthode de Heron;
- Le nombre pi par les différentes méthodes connues (suites qui convergent vers )
- Les intégrales se calculent numériquement par des méthodes analogues;
- La constante d'Euler est définie comme limite d'une suite.
- Les nombres de Bernoulli sont aussi définis comme une suite, bien que pour une fois, ce ne soit pas à sa limite qu'on s'intéresse.
- Même le nombre d'Or peut s'obtenir comme limite d'une suite basée sur les nombres de Fibonacci...
Suites récurrentes
modifierLes suites récurrentes sont celles pour lesquelles dépend de . Mais pour la suite de Fibonacci, la dépendance va plus loin, non seulement le dernier terme intervient, mais également le pénultième.
Suite logistique
modifierLa suite est définie par . Son comportement dépend grandement de la valeur de mais elle est souvent chaotique. Pour calculer ses 20 premiers termes, on peut écrire une simple boucle :
u=0.1
for n in range(20):
u=4*u*(1-u)
print(u)
En effet, une suite récurrente se représente, si n est le temps, par l'affectation d'une variable avec l'image de son ancienne valeur par une fonction (ici, ). Si le premier terme de la suite est une fraction, il en est de même pour tous les termes suivants :
from fractions import *
u=Fraction(1,10)
for n in range(10):
u=4*u*(1-u)
print(u)
Suites arithmétiques et géométriques
modifierUne suite est arithmétique si on passe de chaque terme au suivant en additionnant le même nombre, appelé raison de la suite. L'objet range de Python est donc une suite arithmétique d'entiers.
Par exemple, si on place 2000 € avec des intérêts simples s'élevant à 3 % de 2000, soit 60 € par an, l'évolution du capital pendant 20 ans s'obtient avec ce script :
C=2000
I=2000*3/100
for n in range(1,20):
C+=I
print(C)
Une suite est géométrique si on passe de chaque terme au suivant en multipliant par un même nombre appelé également raison. Par exemple, si on place 2000 € avec des intérêts composés au taux de 2 %, l'évolution du capital année après année est donnée par ce script :
C=2000
for n in range(1,20):
C*=1.02
print(round(C,2))
Le module cmath permet aussi d'étudier les suites géométriques complexes. On constate alors que si le module de la raison est plus petit que 1, la suite tend vers 0, et si le module de la raison est supérieur à 1, la suite tend vers l'infini. C'est bien entendu lorsque le module est égal à 1 qu'il se passe les choses les plus intéressantes...
Suite de Collatz
modifierAlgorithmiquement, la suite de Collatz est intéressante parce que son calcul est basé sur un test de parité, et qu'elle utilise une boucle à condition de sortie :
u=65
while(u>1):
if u%2:
u=3*u+1
else:
u//=2
print(u)
La division par 2 est une division euclidienne, en effet on souhaite que u reste entier (et non flottant) au cours de la boucle.
Suite de Fibonacci
modifierCalcul de la suite
modifierLa récurrence de la suite de Fibonacci est double, avec . Son calcul pose donc un problème algorithmique, puisqu'il faut trois variables (les deux termes à calculer et une variable tampon pour stocker temporairement l'un des deux termes, afin qu'il ne soit pas écrasé par la somme). Ce problème n'existe pas en Python qui permet les affectations simultanées.
a=1
b=1
for n in range(20):
a,b=b,a+b
print(a)
Nombre d'Or
modifierUn problème numériquement intéressant (et c'était la motivation initiale de Fibonacci) est d'étudier le comportement du rapport entre deux termes successifs de la suite de Fibonacci :
a=1
b=1
for n in range(20):
a,b=b,a+b
print(b/a)
print((1+5**0.5)/2)
On constate la convergence vers le nombre d'or.
Suites définies par des sommes
modifierUn exemple
modifierLa suite définie par tend vers 1, il est relativement aisé de le démontrer, et presque aussi facile de le vérifier avec Python:
somme=0
for n in range(1,50):
somme+=1/(n*(n+1))
print(somme)
Un autre exemple
modifierLa suite converge aussi, bien que ce ne soit pas évident en voyant son expression algébrique.
for n in range(1,20):
for k in range(1,n):
somme+=n/(n**2+k)
print(somme)
La constante d'Euler
modifierOn peut la calculer (et vérifier la lenteur de la convergence) avec
from math import *
somme=0
for n in range(1,50):
somme+=1/n
print(somme-log(n))
Calcul de racines carrées
modifierMéthode de Heron
modifierEn constatant que
- Si alors aussi;
- Si alors et vice-versa;
- Par conséquent, on s'attend à ce que la moyenne entre et soit une valeur approchée encore meilleure de ,
on a l'ébauche d'une suite récurrente qui tend vers :
Application
modifieru=1
while(abs(u**2-5)>1e-14):
u=(u+5/u)/2
print(u)
print(5**0.5)
Fonctions en Python
On va traiter un exemple de fonction, issu du Bac STG CGRH Métropole-Réunion de Septembre 2007:
Définition des fonctions
modifierÉnoncé
modifierpour
Et plus bas, dans le même exercice,
pour x appartenant à [5 ; 40].
Fonctions et procédures sous Python
modifierUne fonction est un petit bout de programme Python qui possède un nom (typiquement f), et qui renvoie une valeur (l'image de x, son unique antécédent, par f). Une procédure est une fonction qui ne renvoie rien.
Les fonctions du sujet de Bac ci-dessus peuvent être définies par def en Python:
def C(x):
return x**2+50*x+100.0
def f(x):
return C(x)/x
Si on essaye de calculer l'image de 0 par f, on remarque qu'en lieu et place du calcul, Python affiche un message d'erreur: La fonction f n'est pas definie en 0. On peut aussi la définir directement par
def f(x):
return x+50+100.0/x
On peut aussi définir des fonctions de plusieurs variables (comme le pgcd) ou des fonctions définies par des boucles ou des tests, comme la valeur absolue:
def abs(x):
if x>0:
return x
else:
return -x
Tableau de valeurs
modifierSuite de l'énoncé du Bac 2007
modifierPar la suite, on demande de reproduire et compléter le tableau suivant, arrondi au centième d'euro:
x | 5 | 10 | 20 | 30 | 40 |
f(x) |
Solution en Python
modifierPour remplir le tableau, on peut
- Faire 5 fois un print(f(x)) pour les 5 valeurs de x de l'énoncé;
- Faire une boucle avec les 5 valeurs de x de l'énoncé;
- Faire une boucle avec suffisamment de valeurs de x pour couvrir les 5 valeurs de l'énoncé;
- Transformer f en une fonction qui peut s'appliquer à une liste de valeurs et non à une valeur unique.
La méthode 3 peut se faire avec un itérateur de Python:
for x in range(5,40,5):
print("l\'image de "+str(x)+" par f est "+str(f(x)))
Trop de calculs tuent le calcul! Les images de 25 et 35 par exemple encombrent le tableau. Mais Python permet d'utiliser des itérateurs beaucoup plus souples, et adaptés au problème présent:
for x in [5,10,20,30,40]:
print("l\'image de "+str(x)+" par f est "+str(f(x)))
La méthode 4 peut s'implémenter en Python avec:
print([f(x) for x in [5,10,20,30,40]])
C'est aussi simple que ça!
Arrondis
modifierL'énoncé demandait d'arrondir au centième d'euro près
for x in [5,10,20,30,40]:
print("l\'image de "+str(x)+" par f est "+str(round(f(x),2)))
Représentation graphique
modifierPour représenter graphiquement une fonction comme f ci-dessus, on trace un polygone ayant suffisamment de sommets, et ceux-ci suffisamment proches les uns des autres, pour que la courbe ait l'air courbe. Plusieurs outils permettent de faire du graphisme avec Python mais le plus simple semble être le module tortue.
Avec Turtle
modifierCourbe
modifierLe plus simple est de faire
from turtle import *
setpos(5,f(5))
for x in range(5,40):
setpos(x,f(x))
Difficile de rêver plus simple mais on peut voir trois problèmes:
- Il manque les axes
- Le trait qui relie l'origine au point de coordonnées (5;f(5)) est en trop
- La tortue gêne la visibilité de la figure
Courbe améliorée
modifierPour résoudre le deuxième problème, il suffit de lever le crayon avec penup() avant de commencer le tracé (et de le redescendre pour le tracé lui-même). Pour résoudre le troisième problème, il suffit de rendre la tortue invisible avec hideturtle(). Pour les axes, on peut les tracer avec la tortue:
from turtle import *
penup()
setpos(5,f(5))
pendown()
for x in range(5,40):
goto(x,f(x))
penup()
setpos(0,0)
pendown()
for x in range(0,50,10):
left(90)
forward(2)
backward(4)
forward(2)
right(90)
forward(10)
stamp()
backward(50)
left(90)
for y in range(0,100,10):
left(90)
forward(2)
backward(4)
forward(2)
right(90)
forward(10)
stamp()
Le graphique est correct mais un peu petit et mal cadré.
Avec TkInter
modifierTkInter permet de créer un affichage dans un canevas, lui-même membre d'une fenêtre. Il y a donc plusieurs lignes de Python à écrire avant même de commencer à dessiner. Et l'axe des ordonnées est dirigé vers le bas, ce qui oblige à une transformation des ordonnées. En contrepartie, on a un outil de dessin tout-à-fait correct, et même relativement classique. Et, bien que ce ne soit pas utile dans le cas présent, on peut avoir des boutons, curseurs etc.
On doit donc commencer par
- importer tkInter (s'installer dans l'atelier du peintre)
- créer une fenêtre (monter le chevalet, un cadre...)
- y placer un canevas (tendre la toile sur le chevalet)
- afficher le canevas avec pack (enlever la couverture qui cache le chef-d'œuvre)
- Enfin, dessiner (axes et fonction)
Ce qui peut donner ceci:
from tkinter import *
fenetre=Tk()
graphique=Canvas(fenetre,width=640,height=480)
graphique.pack()
#axe des abscisses
graphique.create_line(20,460,520,460)
for n in range(0,50,10):
graphique.create_line(20+10*n,460,20+10*n,455)
for n in range(0,50,5):
graphique.create_line(20+10*n,460,20+10*n,457)
for n in range(0,50):
graphique.create_line(20+10*n,460,20+10*n,459)
#axe des y
graphique.create_line(20,460,20,60)
for n in range(0,100,10):
graphique.create_line(20,460-4*n,25,460-4*n)
for n in range(0,100,5):
graphique.create_line(20,460-4*n,24,460-4*n)
for n in range(0,100):
graphique.create_line(20,460-4*n,22,460-4*n)
#courbe
for x in range(5,40):
graphique.create_line(20+10*x-10,460-4*f(x-1),20+10*x,460-4*f(x))
Analyse numérique en Python
Fonction
modifierDans ce chapitre, on va effectuer des calculs sur la fonction ; on va appeler cette fonction f. Pour se faciliter la suite, on va créer cette fonction :
def f(x):
return x**2-5
Résolution numérique d'une équation
modifierPour résoudre l'équation f(x)=0, on cherche un intervalle sur lequel on est certain que f s'annule. C'est le cas pour [1;3] parce que f(1) est négatif et f(3) est positif. La méthode de dichotomie vise à resserrer un tel intervalle. On constate ci-dessous que la fonction f est traitée comme une entrée de l'algorithme au même titre que les bornes a et b de l'intervalle :
def zero(f,a,b):
if f(a)*f(b)>0:
print('pas de solution entre '+str(a)+' et '+str(b)+'!')
return 0
while(abs(a-b)>1e-14):
m=(a+b)/2.
if f(m)*f(a)>0:
a=m
else:
b=m
print('la solution de f(x)=0 est '+str(m))
return m
print(zero(f,1,3))
La résolution de l'équation n'est pas terminée, puisque le script ci-dessus n'a donné qu'une seule des deux solutions de cette équation. Par ailleurs, la solution trouvée n'est affichée qu'à près.
Calcul numérique de nombre dérivé
modifierPour calculer le nombre dérivé de f en 5, on va utiliser l'approximation :
def NDer(f,a):
h=1e-10
return (f(a+h)-f(a-h))/(2*h)
print(NDer(f,5))
Calcul numérique d'une intégrale
modifierLa méthode des rectangles dit que où et N est suffisamment grand pour que h soit petit (ci-dessous N=1 000 000) :
def Int(f,a,b):
h=(b-a)/1000000.0
somme=0
for n in range(1000000):
somme+=h*f(a+n*h)
return(somme)
print(Int(f,0,2))
Python et géométrie
Points en Python
L'objet Point est une bonne manière d'aborder la programmation objet. En géométrie repérée, un point est constitué de deux nombres, son abscisse et son ordonnée.
Création de l'objet
modifierLe point de coordonnées (x,y) est, en Python, une classe:
class Point:
def __init__(self,x,y):
self.x=x
self.y=y
Lorsqu'on crée un point, ses coordonnées sont stockées à l'intérieur de l'objet. On note p.x et p.y les coordonnées de p.
Affichage
modifierLa méthode peut ressembler à ceci:
def affichage(self):
return '('+str(self.x)+';'+str(self.y)+')'
mais on peut envisager d'y rajouter des instructions avec TkInter pour réellement dessiner le point sur la figure. Voir à ce sujet le chapitre sur les fonctions.
Avec deux points
modifierLe plus simple quand on a deux points, c'est leur milieu, parce que c'est aussi un point (donc un objet de même nature).
Milieu
modifierLes coordonnées du milieu d'un segment sont les moyennes de celles des extrémités:
def milieu(self,p):
return Point((self.x+p.x)/2,(self.y+p.y)/2)
En se rappelant que l'équivalent en Java de self est this, on remarque une certaine ressemblance avec les codes sources de logiciels de géométrie dynamique:
Tout d'abord, CaRMetal:
setXY((P1.getX() + P2.getX()) / 2, (P1.getY() + P2.getY()) / 2);
Ensuite, GeoGebra:
M.setCoords(
(P.inhomX + Q.inhomX) / 2.0d,
(P.inhomY + Q.inhomY) / 2.0d,
1.0);
Remarque: Ces codes sources sont sous licence GPL ce qui autorise à les citer, au nom de la liberté numéro 1 (celle d'étudier le logiciel) de Richard Stallman.
Vecteur
modifierLe vecteur d'origine A et d'extrémité B, noté , est un vecteur! On le définira donc au chapitre suivant mais il peut servir ici:
def vecteur(self,p):
return Vecteur(p.x-self.x,p.y-self.y)
Distance
modifierPour simplifier l'écriture de la distance AB on va encore utiliser les vecteurs, la distance AB étant égale à :
def distance(self,p):
return self.vecteur(p).norme()
Application au problème
modifierVoici l'objet Point en entier:
from math import *
class Point:
def __init__(self,x,y):
self.x=x
self.y=y
def affichage(self):
return '('+str(self.x)+';'+str(self.y)+')'
def milieu(self,p):
return Point((self.x+p.x)/2,(self.y+p.y)/2)
def vecteur(self,p):
return Vecteur(p.x-self.x,p.y-self.y)
def distance(self,p):
return self.vecteur(p).norme()
Nature de ABC
modifierPour savoir si ABC est isocèle, on peut calculer les longueurs de ses trois côtés:
a=Point(-1,3)
b=Point(5,1)
c=Point(1,5)
print(a.distance(b))
print(a.distance(c))
print(b.distance(c))
Visiblement, ABC n'est pas isocèle. Mais
print(a.distance(b)**2)
print(a.distance(c)**2+b.distance(c)**2)
La réciproque du théorème de Pythagore nous apprend que ABC est rectangle en C, donc d'hypoténuse AB.
Centre du cercle
modifierDonc le cercle circonscrit a pour diamètre [AB], donc pour centre le milieu de [AB]:
m=a.milieu(b)
print(m.affichage())
Rayon du cercle
modifierOn peut donc diviser par 2 la distance AB mais aussi vérifier que M est équidistant de A, B et C:
print(m.distance(a))
print(m.distance(b))
print(m.distance(c))
Figure
modifierOn pourrait utiliser TkInter pour dessiner le tout (y compris le cercle) mais la figure ci-dessous a été faite avec Ruby, ce langage permettant assez facilement de fabriquer un fichier au format svg:
Vecteurs en Python
Géométriquement, un vecteur peut être défini à partir de deux points (son origine et son extrémité) mais aussi par ses coordonnées. C'est le choix qui sera fait ici.
Définition
modifierEncore une fois, on va être classe sur ce coup-là :
class Vecteur:
def __init__(self,x,y):
Lors de son instanciation, un vecteur n'aura donc que deux propriétés : ses coordonnées.
Coordonnées
modifierAbscisse
modifier self.x=x
L'abscisse de u s'obtiendra par u.x
Ordonnée
modifier self.y=y
L'ordonnée de u s'obtiendra par u.y
La classe Vecteur se résume donc pour l'instant à ceci :
class Vecteur:
def __init__(self,x,y):
self.x=x
self.y=y
Affichage
modifierPour afficher un vecteur, on fait comme avec les points (on colle ses coordonnées entre parenthèses, séparées par un point-virgule, après les avoir converties en chaînes de caractères avec str) :
def affichage(self):
return '('+str(self.x)+';'+str(self.y)+')'
Norme
modifierPour calculer la norme d'un vecteur, on utilise le théorème de Pythagore sous la forme de la fonction hypot. Celle-ci doit être importée du module math:
from math import *
La norme du vecteur se calcule donc par
def norme(self):
return hypot(self.x,self.y)
Opérations
modifierPython offre un grand confort, en permettant d'utiliser les signes d'opérations pour les vecteurs. Ainsi, en voyant un signe moins devant un vecteur, Python sait qu'il ne doit pas calculer l'opposé d'un nombre puisque c'est un vecteur, et non un nombre, qui suit ce signe moins. Alors on peut définir (c'est laissé en exercice) une méthode __neg__ pour les vecteurs, et c'est elle qui sera utilisée quand Python verra le signe moins. À moins bien entendu que ce soit une soustraction, auquel cas la méthode __sub__ sera utilisée (cette méthode aussi sera laissée en exercice).
Addition
modifierPour additionner deux vecteurs, on redéfinit la méthode __add__ :
def __add__(self,v):
return Vecteur(self.x+v.x,self.y+v.y)
Alors pour avoir la somme de deux vecteurs u et v, il suffit d'entrer u+v. C'est exactement comme ça que fonctionnent les modules fractions et cmath.
Multiplications
modifierPar un réel
modifierEn multipliant un vecteur par un réel, on obtient un vecteur. Pour ne pas interférer avec la méthode suivante, on va mettre le réel en deuxième et utiliser la méthode multiplication à l'envers qui est notée __rmul__ (comme reverse multiplication) :
def __rmul__(self,r):
return Vecteur(self.x*r,self.y*r)
Pour calculer le triple du vecteur u, on entre 3*u et, pour l'afficher, (3*u).affichage().
Par un autre vecteur
modifierEn multipliant un vecteur par un vecteur, on obtient un nombre, ou scalaire. Aussi cette multiplication est-elle appelée produit scalaire des deux vecteurs (en fait il y a bien une multiplication vectorielle qui donne un vecteur, mais celle-ci est définie pour les vecteurs de l'espace, et ne sera donc pas abordée ici). On peut donc implémenter le produit scalaire de deux vecteurs self et v par
def __mul__(self,v):
return self.x*v.x+self.y*v.y
Pour calculer le produit scalaire de u par v, on entre juste u*v.
Tests
modifierVecteurs colinéaires
modifierPour savoir si deux vecteurs sont colinéaires, on compare deux produits en croix :
def colin(self,v):
return self.x*v.y==self.y*v.x
Ce test retourne un booléen, en entrant u.colin(v).
Vecteurs orthogonaux
modifierPour savoir si deux vecteurs sont orthogonaux, on compare leur produit scalaire à 0 :
def ortho(self,v):
return self*v==0
Exemple
modifierPour l'exercice du chapitre précédent, on peut calculer un produit scalaire pour vérifier que ABC est rectangle :
from math import *
class Vecteur:
def __init__(self,x,y):
self.x=x
self.y=y
def affichage(self):
return '('+str(self.x)+';'+str(self.y)+')'
def norme(self):
return hypot(self.x,self.y)
def __add__(self,v):
return Vecteur(self.x+v.x,self.y+v.y)
def __rmul__(self,r):
return Vecteur(self.x*r,self.y*r)
def __mul__(self,v):
return self.x*v.x+self.y*v.y
def colin(self,v):
return self.x*v.y==self.y*v.x
def ortho(self,v):
return self*v==0
a=Point(-1,3.0)
b=Point(5,1)
c=Point(1,5)
u=c.vecteur(a)
v=c.vecteur(b)
print(u*v)
print(u.ortho(v))
Droites en Python
De façon analogue à la définition par coordonnées des points et vecteurs dans les deux chapitres précédents, on peut définir l'objet droite par une de ses équations. Mais ici, on va définir une droite à partir de deux points. On va donc avoir besoin des classes Point et Vecteur définies dans les deux chapitres précédents. On aura donc aussi besoin d'importer le module math (ou au moins, hypot).
Définition
modifierOn définit donc une classe Droite où l'initialisation place deux points:
class Droite:
def __init__(self,a,b):
self.a=a
self.b=b
Le premier point de la droite d s'obtient par d.a et le second point par d.b.
Vecteurs
modifierVecteur directeur
modifierLe vecteur sera choisi comme vecteur directeur de la droite :
def directeur(self):
return self.a.vecteur(self.b)
On peut s'en servir pour
- Avoir une représentation paramétrique de la droite;
- Faire un test pour savoir si un point est, ou non, sur la droite (par des vecteurs colinéaires).
Ces deux ajouts sont laissés en exercice, le deuxième étant d'ailleurs hors sujet ici puisque c'est une propriété de l'objet Point.
Vecteur normal
modifierLe vecteur normal est choisi de telle façon que son produit scalaire avec le vecteur directeur ci-dessus soit nul:
def normal(self):
return Vecteur(-self.directeur().y,self.directeur().x)
Ce vecteur directeur peut servir à obtenir une équation cartésienne de la droite (AB):
Équations
modifierÉquation cartésienne
modifier def cartesienne(self):
return '('+str(self.normal().x)+')x+('+str(self.normal().y)+')y='+str(self.normal().x*self.a.x+self.normal().y*self.a.y)
On écrit print(d.cartesienne()) pour afficher l'équation cartésienne de d.
Équation réduite
modifierL'équation réduite y=mx+p est acquise dès que le coefficient directeur m et l'ordonnée à l'origine p=y-mx sont acquis.
Coefficient directeur
modifierLe coefficient directeur est le quotient de l'ordonnée du vecteur directeur par son abscisse:
def cd(self):
return self.directeur().y/self.directeur().x
Ordonnée à l'origine
modifier def oalo(self):
return self.a.y-self.cd()*self.a.x
Équation
modifierOn l'obtient à partir des deux nombres précédents:
def reduite(self):
return 'y='+str(self.cd())+'x+('+str(self.oalo())+')'
On l'affiche par print(d.reduite())
Point d'intersection
modifierPour calculer les coordonnées du point d'intersection de deux droites, on doit résoudre un système. Voir à ce sujet le chapitre qui lui est consacré.
Relations entre droites
modifierParallélisme
modifierPour savoir si deux droites sont parallèles, on peut
- Regarder si leurs vecteurs directeurs sont colinéaires (mieux);
- ou comparer leurs coefficients directeurs (moins bien, il faut déjà qu'elles en aient un!)
C'est la première méthode qui sera appliquée ici, bien que la deuxième soit plus rapide:
def parallele(self,d):
return self.directeur().colin(d.directeur())
Orthogonalité
modifierPour savoir si deux droites sont perpendiculaires, on regarde si leurs vecteurs normaux le sont:
def perpendiculaire(self,d):
return self.normal().ortho(d.normal())
Application à l'exercice précédent
modifierEn résumé, l'objet Droite se résume à ceci:
from math import *
class Droite:
def __init__(self,a,b):
self.a=a
self.b=b
def directeur(self):
return self.a.vecteur(self.b)
def normal(self):
return Vecteur(-self.directeur().y,self.directeur().x)
def cartesienne(self):
return '('+str(self.normal().x)+')x+('+str(self.normal().y)+')y='+str(self.normal().x*self.a.x+self.normal().y*self.a.y)
def cd(self):
return self.directeur().y/self.directeur().x
def oalo(self):
return self.a.y-self.cd()*self.a.x
def reduite(self):
return 'y='+str(self.cd())+'x+('+str(self.oalo())+')'
def parallele(self,d):
return self.directeur().colin(d.directeur())
def perpendiculaire(self,d):
return self.normal().ortho(d.normal())
Pour l'exercice des chapitres précédents, on peut vérifier que le triangle ABC est rectangle en C, en regardant si les droites (AC) et (BC) sont perpendiculaires:
a=Point(-1,3.0)
b=Point(5,1)
c=Point(1,5)
d1=Droite(c,a)
d2=Droite(c,b)
print(d1.perpendiculaire(d2))
Résolution de problèmes en Python
Une tortue qui accélère la résolution de problèmes
Depuis la version 2.6, Python possède un module appelé turtle et qui lui permet de faire du graphisme à la LOGO. Outre l'intérêt que peut présenter la consultation de son code source (on y trouve pratiquement tout ce qui est décrit dans les chapitres précédents sur la géométrie), ce module turtle permet de simplifier la résolution de certains problèmes et même d'introduire graphiquement certaines notions mathématiques. Le fil conducteur de ce chapitre est que la tortue LOGO peut mémoriser certaines données de position et se comporte comme une mémoire à la fois plus puissante et moins abstraite que les habituelles variables numériques.
Avant d'utiliser la tortue de Python, on doit l'importer, en faisant
from turtle import *
Ensuite, une connaissance du vocabulaire de situation et de déplacement en anglais peut aider; en voici un échantillon:
- forward ou fd: Pour avancer (l'unité de distance est le pixel)
- backward ou bk: Pour reculer
- left ou lt: Pour tourner à gauche (l'unité d'angle est le degré)
- right ou rt: Pour tourner à droite
- goto pour téléporter la tortue (donner l'abscisse puis l'ordonnée)
- penup() ou pu() pour que les déplacements de la tortue cessent de laisser des traces à l'écran
- pendown() ou pd() pour que les déplacements de la tortue recommencent à laisser des traces
- position() renvoie les coordonnées de la tortue
- home() renvoie la tortue au centre de l'écran (sa position initiale)
- reset() fait pareil mais en effaçant l'écran
- circle dessine un cercle (le rayon est en pixels)
- stamp() donne un coup de tampon sur l'écran, en laissant une empreinte de la tortue
L'écran n'apparaît que lors de l'exécution de la première instruction graphique (un forward par exemple). Sous Windows, il est déconseillé de laisser traîner la souris sur cet écran graphique. Dans les exemples qui suivent, l'export vectoriel au format eps du module TkInter (dont le module turtle hérite) a été utilisé pour produire des figures de meilleure qualité que celles qu'on voit sur l'écran de turtle.
Nombres relatifs
modifierLes nombres (réels) peuvent être représentés par des graduations sur une droite, et donc par les emplacements de la tortue à l'écran.
Nombres positifs
modifierAddition
modifierPour représenter l'addition de 21 et 34, on peut tout simplement entrer
from turtle import *
forward(21)
forward(34)
print(position())
Ce qui oblige à ignorer une information superflue (l'ordonnée de la tortue). La variante suivante permet d'éviter cela:
reset()
forward(21)
forward(34)
print(distance(0,0))
Soustraction
modifierPour soustraire 21 à 34, il suffit de faire reculer la tortue au lieu de la faire avancer:
reset()
forward(34)
backward(21)
print(position())
Si on intervertit l'amplitude des mouvements, on découvre que Python choisit d'afficher négativement une position à gauche de l'origine:
reset()
forward(21)
backward(34)
print(position())
Assez naturellement, on est amené à poser 21-34=-13: Découverte expérimentale des nombres négatifs...
Nombres négatifs
modifierUne fois qu'on a vu des nombres négatifs, on peut chercher comment réaliser des opérations dessus:
Addition
modifierPour additionner deux nombres négatifs, on peut faire
reset()
backward(34)
backward(21)
position()
Tout ceci permet assez rapidement d'explorer les différents autres cas de figure (deux cas différents pour la somme de deux nombres de signes différents). Puis la découverte spontanée du fait que les deux instructions suivantes ont le même effet:
forward(-34)
backward(34)
Ce qui facilite grandement l'exploration de la soustraction de deux nombres relatifs:
Soustraction
modifierPour calculer 34-(-21), on peut faire
reset()
forward(34)
backward(-21)
position()
Pour l'apprentissage des opérations sur les nombres négatifs, turtle constitue un outil expérimental intéressant à explorer.
Angles orientés
modifierDe même, les deux instructions suivantes ont le même effet (rotation de 60° vers la gauche):
left(60)
right(-60)
mais ce n'est nullement évident pour des lycéens qui n'ont jamais fait ce genre de manipulation, surtout depuis que la notion de rotation a totalement disparu de l'enseignement des mathématiques. Pourtant le module turtle permet de visualiser l'addition des angles et d'introduire des notions comme celle d'angles complémentaires ou supplémentaires avec
left(30)
left(60)
Fonctions
modifierOn a vu dans un chapitre précédent comment la module turtle permet de représenter graphiquement une fonction.
Statistique
modifierChute d'une bille sur la planche de Galton
modifierLa planche de Galton réalise une marche aléatoire de dimension 1 (le mouvement vertical de la bille n'ayant aucune influence sur le numéro de la case où elle aboutit). On peut donc simuler le mouvement d'une bille avec ce script:
from turtle import *
from random import *
for n in range(24):
if random()<0.5:
forward(1)
else:
backward(1)
print(position())
On peut améliorer ce script en utilisant randrange qui va de -1 à 1 (donc 2 exclu) par pas de 2, ce qui économise un test:
from turtle import *
from random import *
for h in range(24):
forward(randrange(-1,2,2))
print(position())
Statistiques sur 100 billes
modifierPour effectuer des statistiques sur 100 billes, on a intérêt à accélérer la tortue, avec
from turtle import *
from random import *
speed=1000
hideturtle()
penup()
Ensuite on crée un tableau d'effectifs pour simuler le bas de la planche de Galton:
effectifs=[0 for x in range(-24,25)]
Après ça il n'y a plus qu'à remplir le tableau en recommençant 100 fois l'expérience précédente (lancer d'une bille):
from turtle import *
from random import *
speed=0
hideturtle()
penup()
for n in range(100):
home()
for h in range(24):
forward(randrange(-1,2,2))
effectifs[int(xcor()]+=1
Ce script, bien qu'assez court, met du temps à s'exécuter (de l'ordre d'une minute). Pour l'accélérer, on peut ajouter un degré d'abstraction en n'utilisant pas la tortue. En effet, chaque pas est égal à -1 ou 1 au hasard, donc d'après ce qu'on a vu au début de ce chapitre (opérations sur les nombres relatifs), on ne fait qu'additionner 24 nombres égaux à 1 ou -1, ce qui donne ce script:
from random import *
effectifs=[0 for x in range(-24,25)]
for n in range(100):
effectifs[sum(randrange(-1,2,2) for h in range(24))]+=1
Cette fois-ci, l'effet est presque instantané.
Dessin de l'histogramme
modifierUne fois le tableau d'effectifs rempli, le module turtle peut le représenter graphiquement sous forme d'un polygone des effectifs. Comme la méthode précédente est très rapide et que les effectifs des nombres impairs sont nuls, on va plutôt utiliser la variante suivante, avec 256 cases et 1000 essais:
from turtle import *
from random import *
effectifs=[0 for x in range(256)]
for n in range(1000):
effectifs[sum(randrange(2) for h in range(256))]+=1
reset()
for x in range(256):
goto(x,effectifs[x])
On obtient alors un histogramme de ce genre (la tortue est encore visible à droite):
Fractales
modifierLa courbe de Von Koch est classiquement définie par la récursivité. Mais elle n'est pas nécessaire si on utilise une expression régulière. En fait, on peut écrire un script Python qui produit un script Python, puis exécuter celui-ci!
Triangle de départ
modifierComme le script qui va dessiner le triangle fractal sera assez long, on va utiliser des abréviations: fd au lieu de forward, lt au lieu de left et rt au lieu de right. Alors pour dessiner un triangle on peut faire ceci:
from turtle import *
fd(100); rt(120); fd(100); rt(120); fd(100); rt(120)
Ou mieux, en stockant ce programme en Python dans une variable programme:
from turtle import *
programme='fd(100); rt(120); fd(100); rt(120); fd(100); rt(120)'
exec(programme)
Dans un premier temps, on va abréger encore plus, en notant chaque instruction de ce programme par une seule lettre:
- A (comme avance) pour fd(100);
- p (comme plus) pour lt(60);
- m (comme moins) pour rt(120).
La traduction se fait par une RegExp, qui, tel un chien de douane, cherche toutes les occurences d'une lettre, et les remplace par le texte correspondant.
Alors le programme pour créer un programme qui dessine un triangle devient:
from turtle import *
from re import *
programme='AmAmAm'
programme=sub('A','fd(100); ',programme)
programme=sub('m','rt(120); ',programme)
exec(programme)
Le remplacement des lettres mnémotechniques par des instructions en Python est à l'image de ce que fait un compilateur comme celui de Python. Avec ça, au moins, la recette pour dessiner un triangle est facile à retenir: avancer; tourner; avancer; tourner; avancer; tourner, étant entendu que chaque fois qu'on avance, c'est de 100 pixels, et chaque fois qu'on tourne, c'est de 120° vers la droite.
Modification du script
modifierPour transformer le triangle en flocon, on doit remplacer chaque instruction avancer par la séquence avancer; gauche; avancer; droite; avancer; gauche; avancer. Du moment que chaque fois qu'on avance, c'est du même nombre de pixels (par exemple 81) et chaque fois qu'on tourne à gauche, c'est de 60° et chaque fois qu'on tourne à droite, c'est de 120°. Pour obtenir cet effet, il suffit de remplacer chaque A par ApAmApA:
from turtle import *
from re import *
programme='AmAmAm'
programme=sub('A','ApAmApA',programme)
programme=sub('A','fd(81); ',programme)
programme=sub('m','rt(120); ',programme)
programme=sub('p','lt(60); ',programme)
exec(programme)
Ce script dessine bien une étoile:
Dessin du triangle de Von Koch
modifierPour finir le dessin du flocon fractal, il suffit d'itérer le remplacement de chaque A par ApAmApA:
from turtle import *
from re import *
programme='AmAmAm'
for n in range(4):
programme=sub('A','ApAmApA',programme)
programme=sub('A','fd(2); ',programme)
programme=sub('m','rt(120); ',programme)
programme=sub('p','lt(60); ',programme)
exec(programme)
Ce script dessine ceci en 9 lignes de Python (la tortue donne une idée de l'échelle):
Voir aussi
modifierSur ce site, il y a un wikibook sur LOGO dont beaucoup d'idées sont transposables ici.
Résolution de systèmes en Python
Bien entendu, on va noter x le prix d'un xylophone, et y le prix d'un youkoulélé. Les données de l'énoncé se traduisant algébriquement par 2y=3x+1 et x+y=8.
On va donc voir comment résoudre le système de deux équations à deux inconnues suivant :
Méthode graphique
modifierUne méthode simple (surtout ici puisque la solution est formée d'entiers) consiste à tracer les deux droites d'équations respectives et et de lire sur le graphique les coordonnées de leur point d'intersection. Python tout seul ne sait pas faire ça mais PyKig lance le logiciel Kig et le laisse faire le travail, à condition de lui fournir les éléments nécessaires pour faire les constructions géométriques : Les coordonnées de points de chaque droite. Or la droite d'équation passe par les points de coordonnées et (intersections avec les axes de coordonnées) qui sont constructibles par Kig. Le script suivant colorie en bleu la première droite (le milieu M sert de point d'ancrage pour son équation réduite) et en rouge la seconde droite, puis affiche en mauve leur point d'intersection :
e1=[3.0,-2.0,-1.0]
e2=[1.0,1.0,8.0]
A=Point(e1[2]/e1[0],0,HIDDEN)
B=Point(0,e1[2]/e1[1],HIDDEN)
d1=Line(A,B)
d1.setcolor("blue")
M=MidPoints(A,B,HIDDEN)
t=Text(M,Equation(d1),0)
t.setcolor("blue")
C=Point(e2[2]/e2[0],0,HIDDEN)
D=Point(0,e2[2]/e2[1],HIDDEN)
d2=Line(C,D)
d2.setcolor("red")
N=MidPoints(C,D,HIDDEN)
u=Text(N,Equation(d2),0)
u.setcolor("red")
I=LineLineIntersection(d1,d2)
I.setcolor("magenta")
En enregistrant son contenu dans un fichier appelé system.kpy, et en exécutant dans une console ceci :
pykig.py system.kpy
On obtient, après lancement inopiné de Kig (à condition bien sûr que celui-ci soit installé, ce qui est très facile sous Linux, beaucoup moins sous les autres systèmes), la figure suivante (où il a été nécessaire de créer un Label donnant les coordonnées du point d'intersection, ce que PyKig ne sait pas faire) :
On lit les coordonnées du point d'intersection, qui constituent la solution du système. Il est possible de déplacer à la souris les points A, B, C et D et voir les coordonnées du point d'intersection se mettre à jour en temps réel : En effet Kig est un logiciel de géométrie dynamique.
Cette méthode n'est pas applicable à des systèmes de plus de deux inconnues mais elle est très visuelle pour des systèmes tels que pour lesquels on voit non seulement qu'il n'y a pas de solution, mais également pourquoi il n'y a pas de solution.
Si l'une des droites est verticale ou horizontale, le fichier ci-dessus doit être modifié avec des tests gérant ces cas. Ceci est laissé en exercice.
Méthode itérative
modifierSi on sait que x et y sont entiers naturels, on peut résoudre le système par une double boucle sur x et y :
solutions=[(x,y) for x in range(100) for y in range(200) if 3*x-2*y==-1 and x+y==8]
print(solutions)
Pour "x"et"y" appartenant à R je propose le script suivant :
def systeme (a1,b1,c1,a2,b2,c2):
x=float() y=float() 'a1*x + b1*y =c1,a2*x + b2*y =c2' if a1*b2-a2*b1==0: print('Pas de solution') else: y=(c2*a1-c1*a2)/(a1*b2-a2*b1) x=(c1-b1*y)/a1 print('x =',round(x,2),"",'y =',round(y,2))
Méthode de Cramer
modifierLa méthode de Cramer est également implémentable en Python :
def affiche(e):
print(str(e[0])+'x+('+str(e[1])+')y='+str(e[2]))
def resoudre(e1,e2):
determinant=e1[0]*e2[1]-e1[1]*e2[0]
if determinant==0:
print('Pas de solution unique')
else:
x=(e1[2]*e2[1]-e1[1]*e2[2])/determinant
y=(e1[0]*e2[2]-e1[2]*e2[0])/determinant
print('La solution est ('+str(x)+','+str(y)+')')
e=[3,-2,-1]
f=[1,1,8]
affiche(e)
affiche(f)
resoudre(e,f)
Mais le module fractions permet d'avoir la valeur exacte de la solution chaque fois que les coefficients du système sont entiers :
from fractions import *
def resoudre(e1,e2):
determinant=e1[0]*e2[1]-e1[1]*e2[0]
if determinant==0:
print('Pas de solution unique')
else:
x=Fraction(e1[2]*e2[1]-e1[1]*e2[2],determinant)
y=Fraction(e1[0]*e2[2]-e1[2]*e2[0],determinant)
print('La solution est ('+str(x)+','+str(y)+')')
e=[3,-2,-1]
f=[1,1,8]
resoudre(e,f)
Avec NumPy
modifierLe module NumPy permet de faire du calcul avec Python, et même en étant à la fois rapide et précis (parce que ce module est précompilé).
Installation
modifierLes principales distributions Linux sont fournies avec. Sur Windows, on peut le télécharger sur http://sourceforge.net/projects/numpy/files/.
Une fois le .zip décompressé, il faut l'installer en tapant dans la console :
python setup.py install
En cas d'erreur :
- ImportError “No Module named Setuptools” : sauvegarder le fichier https://bootstrap.pypa.io/ez_setup.py dans le dossier Python, puis l'installer :
python ez_setup.py install
.- Microsoft Visual C++ 9.0 is required, le télécharger sur https://www.microsoft.com/en-us/download/details.aspx?id=44266
Utilisation
modifierPour résoudre le système, entrer le script suivant :
from numpy import *
A=matrix([[3,-2],[1,1]])
B=matrix([[-1],[8]])
solution=linalg.solve(A,B)
print(solution)
Triplets pythagoriciens en Python
L'énoncé du problème est simple, sa solution avec Python aussi:
Trouver tous les triplets pythagoriciens (x,y,z) tels que ; autrement dit, on demande les triangles rectangles de périmètre 1000 dont les côtés sont entiers.
En considérant x, y et z comme classés dans l'ordre croissant, on va faire une boucle sur y (le plus grand des côtés de l'angle droit), et à l'intérieur de celle-ci, une autre boucle sur x. Enfin on calcule l'hypoténuse z, puis le périmètre du triangle. Et si celui-ci vaut 1000, on affiche le triplet:
from math import hypot
for y in range(1000):
for x in range(y):
z=hypot(x,y)
if x+y+z==1000:
print(x,y,z)
Le script nous apprend qu'il n'y a qu'un seul triplet pythagoricien de somme 1000.
Systèmes congruentiels en Python
Rallye mathématique de la Réunion
modifierSujet 2005, Exercice 2
modifierOn demande combien d'élèves il y a au total.
On peut utiliser une boucle avec plein de tests pour détecter les solutions mais aussi les ensembles. C'est ce qu'on va faire ici.
Construction de l'ensemble 2 par 2
modifierOn commence par construire la liste des nombres entre 2800 et 2900 tels que 2 par 2, il en reste 1: Ces nombres constituent une suite arithmétique de raison 2.
s2=[n for n in range(2800,2900) if n%2==1]
print(len(s2))
À ce stade, il y a beaucoup de candidats possibles:
{2817, 2819, 2821, 2823, 2825, 2827, 2829, 2831, 2833, 2835, 2837, 2839, 2841, 2843, 2845, 2847, 2849, 2851, 2853, 2855, 2857, 2859, 2861, 2863, 2865, 2867, 2869, 2871, 2873, 2875, 2877, 2879, 2881, 2883, 2885, 2887, 2889, 2891, 2893, 2895, 2897, 2899, 2801, 2803, 2805, 2807, 2809, 2811, 2813, 2815}
Épuration 3 par 3
modifierMaintenant on va construire un ensemble analogue s3 pour les nombres tels que, 3 par 3, il en reste 1, choisis parmi les nombres de s2:
s3=[n for n in s2 if n%3==1]
print(s3)
Il en reste déjà moins:
{2881, 2851, 2821, 2887, 2857, 2827, 2893, 2863, 2833, 2899, 2869, 2803, 2839, 2809, 2875, 2845, 2815}
Cas du nombre 4
modifierCette fois-ci, dans s4 on met les nombres précédents tels que pris 4 par 4, il en reste 1:
s4=[n for n in s3 if n%4==1]
print(s4)
{2881, 2821, 2857, 2893, 2833, 2869, 2809, 2845}
Ça se précise!
Avec 5
modifierMaintenant on va de 5 en 5:
s5=[n for n in s4 if n%5==1]
print(s5)
{2881, 2821}
Il ne reste que deux nombres à tester!
Avec 6
modifierOn pourrait les tester l'un après l'autre, mais aussi faire pareil qu'avant (on y prend goût !):
s6=[n for n in s5 if n%6==1]
print(s6==s5)
Aucun changement, toujours deux solutions:
{2881, 2821}
Dernière étape
modifierAutant continuer sur la même voie:
solutions=[n for n in s6 if n%7==0]
print(solutions)
Ce qui donne la réponse (on constate qu'elle est unique) à la question de l'énoncé.
Sujet 2007, Exercice 3
modifierCe soir, veille de marché, il est perplexe.
• S’il met ses œufs dans des emballages de 6, il en reste 2.
• S’il utilise des emballages de 10, il en reste encore 2.
• Il me faudrait, dit-il, des emballages de 8 pour tout contenir exactement.
On demande combien il y a d'œufs en tout.
Par 6
modifierOn commence par examiner les nombres tels que, pris 6 par 6, il en reste toujours 2:
s6=[n for n in range(40,200) if n%6==2]
print(len(s6))
26 nombres à examiner, courage!
Par 10
modifiers10=[n for n in s6 if n%10==2]
print(s10)
Plus que 5 nombres à examiner!
Par 8
modifiersolutions=[n for n in s10 if n%8==0]
print(solutions)
Comme une seule valeur de n est affichée, c'est la seule qui a réussi tous les tests, et c'est donc l'unique solution au problème.
Freudenthal en Python
Le problème de Freudenthal est intéressant à traiter en Python parce qu'il peut se résoudre en manipulant des objets (tableaux et ensembles). Le problème est d'ailleurs intéressant en soi parce qu'il porte sur la logique épistémique.
- Polly : je ne sais pas qui sont x et y.
- Sam : je savais que tu ne savais pas !
- Polly : alors je sais qui sont x et y.
- Sam : alors je sais aussi !
Le but du problème est donc d'utiliser les connaissances qu'on a sur les connaissances de Polly et Sam pour savoir si on peut savoir quels sont x et y.
Un Outil
modifierPour se faciliter la suite, on va créer une fonction Python qui, à chaque somme n, associe la liste des produits qui ont donné cette somme :
def prod(n):
p=[]
for k in range(2,n-2):
p.append(k*(n-k))
return p
Alors prod(8) donne la liste [12, 15, 16, 15, 12] parce que 8 peut s'écrire 2+6, 3+5, 4+4, 5+3 ou 6+2 et que les produits correspondants sont 12, 15 et 16 (certains apparaissant deux fois).
Première affirmation
modifierL'affirmation apporte une information : Le produit donné à Polly peut s'obtenir de plusieurs manières, sinon Polly connaîtrait les facteurs. Pour exploiter cette information, on va commencer par fabriquer l'énorme liste des produits possibles, puis ne garder que ceux qui apparaissent au moins deux fois dans la liste :
produits=[x*y for y in range(3,100) for x in range(2,y) if x+y<=100]
print(len(produits))
polly=[p for p in produits if produits.count(p)>=2]
print(len(polly))
Ceci dit, il en reste encore pas mal...
Deuxième affirmation
modifierSi Sam sait que Polly ne sait pas, c'est parce que quelle que soit la décomposition en somme d'entiers de celui qu'on lui a dicté, le produit correspondant est dans la liste précédente. Sinon Polly aurait pu l'entendre et aurait alors su quels en sont les facteurs. Sam ne va donc garder que les sommes n pour lesquelles la liste prod(n) calculée avec la fonction ci-dessus ne contient que des éléments de la liste polly, donc si leur intersection a une longueur maximale (soit n-4) :
sam=[n for n in range(4,100) if len([p for p in prod(n) if p in polly])==n-4]
print(sam)
Mais le nombre 4 qui apparaît est une erreur, dûe à ce que la liste est vide, donc de longueur 4-4=0. Or 4 ne peut s'écrire que 2+2 et x et y sont supposés différents. Donc on va plutôt commencer par 5 :
sam=[n for n in range(5,100) if len([p for p in prod(n) if p in polly])==n-4]
print(sam)
On voit alors apparaître la liste des sommes que Sam a pu somme toute ouïr :
[11, 17, 23, 27, 29, 35, 37, 41, 47, 53]
Troisième affirmation
modifierLa dernière affirmation de Sam lève toute ambigüité chez Polly. Mais quel genre d'ambiguïté peut-il y avoir encore? Par exemple, un des produits associés à 11 est 30 (car 11=6+5) et 30 est aussi un produit associé à 17 (car 17=15+2). Si le produit 30 avait été confié à Polly, l'ambiguïté en question resterait présente. Polly va donc enlever aux listes prod(11), prod(17) etc. les doublons. On va donc créer la liste des doublons avec
doublons=[]
for p in sam:
for q in sam:
if q<>p:
doublons+=([r for r in prod(p) if r in prod(q)])
Puis on va enlever à chaque liste des produits (de 11, de 17 etc.) chaque doublon. Si Polly connaît la somme de Sam, c'est parce que dans la liste sam, il ne restera qu'une liste de produits contenant exactement deux produits ( et ). Il ne reste alors plus qu'à la chercher pour en savoir autant que Polly :
solutions=[p for p in sam if len([r for r in prod(p) if r not in doublons])==2]
print(solutions)
On connaît donc la somme de Sam.
Quatrième affirmation
modifierPuisque la somme de Sam vaut 17, on recommence l'étape précédente avec 17 seulement : Chercher la liste des produits de 17, doublons enlevés :
print([r for r in prod(17) if r not in doublons])
On connaît maintenant le produit de Polly.
Recherche de x et y
modifierTrouver deux nombres dont le produit est 52 et la somme 17 est un problème classique sur les équations du second degré. Mais Python permet aussi de les trouver avec une double boucle (puisqu'on sait que ces nombres sont entiers) :
print([(x,y) for y in range(100) for x in range(y) if x+y==17 and x*y==52])
La solution (x,y) apparaît alors entre crochets.
GFDL | Vous avez la permission de copier, distribuer et/ou modifier ce document selon les termes de la licence de documentation libre GNU, version 1.2 ou plus récente publiée par la Free Software Foundation ; sans sections inaltérables, sans texte de première page de couverture et sans texte de dernière page de couverture. |