Aller au contenu

Mathématiques avec Python et Ruby/Version imprimable1

Un livre de Wikilivres.

Ceci est la version imprimable de Mathématiques avec Python et Ruby.
  • Si vous imprimez cette page, choisissez « Aperçu avant impression » dans votre navigateur, ou cliquez sur le lien Version imprimable dans la boîte à outils, vous verrez cette page sans ce message, ni éléments de navigation sur la gauche ou en haut.
  • Cliquez sur Rafraîchir cette page pour obtenir la dernière version du wikilivre.
  • Pour plus d'informations sur les version imprimables, y compris la manière d'obtenir une version PDF, vous pouvez lire l'article Versions imprimables.


Mathématiques avec Python et Ruby

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

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. Une copie de cette licence est incluse dans l'annexe nommée « Licence de documentation libre GNU ».

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:

  1. int pour les nombres entiers;
  2. long pour les grands entiers;
  3. fraction pour les fractions (quotients d'entiers par des entiers);
  4. decimal et float pour les réels (en réalité des nombres décimaux);
  5. 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

[modifier | modifier le wikicode]

Ainsi, 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())

Addition, soustraction et multiplication

[modifier | modifier le wikicode]

Les 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)

On 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())

Il 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):

Dans 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

[modifier | modifier le wikicode]

Le 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.


Deux 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.

Beaucoup 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

[modifier | modifier le wikicode]

En Python comme en algèbre, on effectue dans l'ordre

  1. Les parenthèses
  2. Les fonctions (comme l'élévation à une puissance)
  3. Les multiplications et divisions
  4. 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.

Un 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 :

  1. 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);
  2. Lorsqu'on dit qu'il est 8 heures 13, c'est que depuis minuit, il est passé exactement heures (soit 493 minutes).
  3. 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...
  4. 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".
  5. 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

[modifier | modifier le wikicode]

Pour 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 ...

Les 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

[modifier | modifier le wikicode]

L'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)

Pour obtenir l'inverse d'une fraction, on divise 1 par celle-ci :

from fractions import *
a=Fraction(5,4)
print(1/a)

La somme de deux fractions est une fraction :

from fractions import *
a=Fraction(34,21)
b=Fraction(21,13)
print(a+b)

La différence de deux fractions est une fraction :

from fractions import *
a=Fraction(34,21)
b=Fraction(21,13)
print(a-b)


Multiplication

[modifier | modifier le wikicode]

Le produit de deux fractions est une fraction :

from fractions import *
a=Fraction(34,21)
b=Fraction(21,13)
print(a*b)


Le 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)

Si 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)


Réduite de Farey

[modifier | modifier le wikicode]

La 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

[modifier | modifier le wikicode]

Une 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

  1. Il a fallu importer le module math qui contient cette fonction ceil;
  2. 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

[modifier | modifier le wikicode]

Si 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

[modifier | modifier le wikicode]

Un 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.

Une 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))

Pour "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

[modifier | modifier le wikicode]

Le 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

[modifier | modifier le wikicode]

Deux 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

[modifier | modifier le wikicode]

Les 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

[modifier | modifier le wikicode]

Puissances et exposants

[modifier | modifier le wikicode]

Comme 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.

Le 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

[modifier | modifier le wikicode]

Pour 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

[modifier | modifier le wikicode]

On 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

[modifier | modifier le wikicode]

Dans 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)


Les 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

[modifier | modifier le wikicode]

Les 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

[modifier | modifier le wikicode]

Pour 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))

Avec cmath, on peut appliquer certaines fonctions de la variable réelle à des complexes.

Exponentielles

[modifier | modifier le wikicode]

Pour 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))

On 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

[modifier | modifier le wikicode]

Le 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))


Le 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

On 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

Toutes 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.

Pour 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).

	def __neg__(self):
		return Quaternion(-self.a,-self.b)

En écrivant -q, on aura désormais l'opposé de q.

	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.

	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é.

Pour 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)
	def __sub__(self,other):
		return Quaternion(self.a-other.a,self.b-other.b)

Multiplication

[modifier | modifier le wikicode]

Le 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)
Multiplication par un réel
[modifier | modifier le wikicode]

Pour 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.

Le 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 : .

Mê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 q et c sont des quaternions, mène à des ensembles de Mandelbrot quaternionniques ([1])

Voici 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.

Ce 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

[modifier | modifier le wikicode]
from math import hypot

class Octonion:
	
	def __init__(self,a,b):
		self.a=a
		self.b=b

Comme 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!

Les fonctions sur les octonions se définissent presque comme celles sur les quaternions, Cayley-Dickson oblige :

	def __neg__(self):
		return Octonion(-self.a,-self.b)
	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.

	def conjugate(self):
		return Octonion(self.a.conjugate(),-self.b)
	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.

	def __sub__(self,other):
		return Octonion(self.a-other.a,self.b-other.b)

Multiplication

[modifier | modifier le wikicode]
	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))

Grâ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 | modifier le wikicode]
	def __rmul__(self,k):
		return Octonion(k*self.a,k*self.b)
Quotient de deux octonions
[modifier | modifier le wikicode]
	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.

Comme 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.


La 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 *
  • 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 | modifier le wikicode]

Évènements certain et impossible

[modifier | modifier le wikicode]

L'é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.

On 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

[modifier | modifier le wikicode]

Cette 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 | modifier le wikicode]

Évènements simultanés

[modifier | modifier le wikicode]

L'é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.

univers={1,2,3,4,5,6}
petit={1,2,3,4}
pair={2,4,6}
print(petit&pair)


Avec les cartes

[modifier | modifier le wikicode]
print(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.

De même l'évènement "A ou B" se note , et en Python, le symbole pipe (trait vertical). Python enlève automatiquement les doublons.

univers={1,2,3,4,5,6}
petit={1,2,3,4}
pair={2,4,6}
print(petit|pair)

Avec les cartes

[modifier | modifier le wikicode]
print(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

[modifier | modifier le wikicode]

Pour calculer le contraire d'un évènement, on le soustrait à l'univers.

univers={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

[modifier | modifier le wikicode]
print(univers-figure)
print(univers-pique)


On 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.

def 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

[modifier | modifier le wikicode]
def 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.

Versions jusqu'à la 2.7

[modifier | modifier le wikicode]

Python 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:

  1. Vérifier que les modules 30269, 30307 et 30323 sont premiers.
  2. Vérifier que les multiplicateurs 171, 172 et 170 sont primitifs modulo leurs modules respectifs.
  3. 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 ).
  4. 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...

Python 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.

Toutes ces merveilles se trouvant dans le module random, on doit charger celui-ci pour faire des simulations avec Python:

from random import *


Variables uniformes

[modifier | modifier le wikicode]

Le 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))


Pour 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

[modifier | modifier le wikicode]

La 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

[modifier | modifier le wikicode]

On 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

[modifier | modifier le wikicode]

En 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.

Il le définit ainsi:

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),

Euler montre que

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

[modifier | modifier le wikicode]

Un 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

[modifier | modifier le wikicode]

Pour 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)

Pour é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

[modifier | modifier le wikicode]

On 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)

On tire 5 cartes d'un jeu de 32. Quelle est la probabilité des évènements suivants:

  1. On a une couleur (les 5 cartes sont de la même couleur);
  2. 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

[modifier | modifier le wikicode]

On 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

[modifier | modifier le wikicode]

Pour 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,

  1. On se fixe des probabilités a priori de se tromper dans ses estimations ou son test;
  2. On cherche (difficile) pour quels choix des paramètres (intervalle ou estimateur) ces probabilités sont atteintes;
  3. 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.

Pour 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))

Pour 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.

La 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.

L'écart-type est défini comme la racine carrée de la variance:

def ecartype(tableau):
    return variance(tableau)**0.5

print(ecartype(donnees))

On 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

[modifier | modifier le wikicode]

Pour des fréquences

[modifier | modifier le wikicode]
Une situation fictive
Sur les 100 000 électeurs d'une ville, 43 000 s'apprêtent à voter pour le maire sortant, mais celui-ci ne le sait pas. Alors il commande un sondage, et l'institut de sondage constitue un échantillon de 100 habitants sur lesquels 52 disent vouloir voter pour le maire. Fou de joie, celui-ci achète le champagne pendant la campagne, avant les élections, et se fait battre lors des élections. Il accuse l'institut de sondage d'avoir triché. Qu'en penser?


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

[modifier | modifier le wikicode]

Là 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

[modifier | modifier le wikicode]

En 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

[modifier | modifier le wikicode]

Euler 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:

  1. Pour calculer un nombre , on invente une suite telle que ;
  2. On estime que par exemple 1000 est suffisamment proche de pour que .
  3. On calcule alors (en général par récurrence).
  4. On considère donc qu'on a fini de calculer .

C'est ainsi par exemple, qu'on calcule

  1. des racines carrées par la méthode de Heron;
  2. Le nombre pi par les différentes méthodes connues (suites qui convergent vers )
  3. Les intégrales se calculent numériquement par des méthodes analogues;
  4. La constante d'Euler est définie comme limite d'une suite.
  5. 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.
  6. Même le nombre d'Or peut s'obtenir comme limite d'une suite basée sur les nombres de Fibonacci...

Suites récurrentes

[modifier | modifier le wikicode]

Les 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

[modifier | modifier le wikicode]

La 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

[modifier | modifier le wikicode]

Une 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

[modifier | modifier le wikicode]

Algorithmiquement, 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

[modifier | modifier le wikicode]

Calcul de la suite

[modifier | modifier le wikicode]

La 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)

Un 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

[modifier | modifier le wikicode]

La 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

[modifier | modifier le wikicode]

La 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

[modifier | modifier le wikicode]

On 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

[modifier | modifier le wikicode]

Méthode de Heron

[modifier | modifier le wikicode]

En constatant que

  1. Si alors aussi;
  2. Si alors et vice-versa;
  3. 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 :

u=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 | modifier le wikicode]
Extrait du sujet de l'exercice 3:
Une entreprise produit des appareils électroménagers. Le coût horaire de production de x appareils est donné en euros par :

pour

Et plus bas, dans le même exercice,

le coût moyen est défini par
Le coût moyen de production d’un objet est égal à

pour x appartenant à [5 ; 40].

Fonctions et procédures sous Python

[modifier | modifier le wikicode]

Une 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

[modifier | modifier le wikicode]

Suite de l'énoncé du Bac 2007

[modifier | modifier le wikicode]

Par 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

[modifier | modifier le wikicode]

Pour remplir le tableau, on peut

  1. Faire 5 fois un print(f(x)) pour les 5 valeurs de x de l'énoncé;
  2. Faire une boucle avec les 5 valeurs de x de l'énoncé;
  3. Faire une boucle avec suffisamment de valeurs de x pour couvrir les 5 valeurs de l'énoncé;
  4. 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!

L'é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

[modifier | modifier le wikicode]

Pour 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.

Le 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:

  1. Il manque les axes
  2. Le trait qui relie l'origine au point de coordonnées (5;f(5)) est en trop
  3. La tortue gêne la visibilité de la figure

Courbe améliorée

[modifier | modifier le wikicode]

Pour 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é.

TkInter 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

  1. importer tkInter (s'installer dans l'atelier du peintre)
  2. créer une fenêtre (monter le chevalet, un cadre...)
  3. y placer un canevas (tendre la toile sur le chevalet)
  4. afficher le canevas avec pack (enlever la couverture qui cache le chef-d'œuvre)
  5. 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

Dans 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

[modifier | modifier le wikicode]

Pour 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é

[modifier | modifier le wikicode]

Pour 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

[modifier | modifier le wikicode]

La méthode des rectangles dit que 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.


Voici l'énoncé de l'exercice:
Dans un repère orthonormé, on considère , et . Calculer les distances AB, AC et BC et en déduire la nature du triangle ABC. Puis en déduire les coordonnées du centre de son cercle circonscrit.

Création de l'objet

[modifier | modifier le wikicode]

Le 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.

La 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

[modifier | modifier le wikicode]

Le 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).

Les 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.

Le 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)

Pour 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

[modifier | modifier le wikicode]

Voici 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()

Pour 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

[modifier | modifier le wikicode]

Donc 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

[modifier | modifier le wikicode]

On 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))

On 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.

Encore 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.

        self.x=x

L'abscisse de u s'obtiendra par u.x

        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

Pour 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)+')'

Pour 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)

Python 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).

Pour 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

[modifier | modifier le wikicode]

En 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

[modifier | modifier le wikicode]

En 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.

Vecteurs colinéaires

[modifier | modifier le wikicode]

Pour 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

[modifier | modifier le wikicode]

Pour savoir si deux vecteurs sont orthogonaux, on compare leur produit scalaire à 0 :

    def ortho(self,v):
        return self*v==0

Pour 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).

On 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.

Vecteur directeur

[modifier | modifier le wikicode]

Le vecteur sera choisi comme vecteur directeur de la droite :

    def directeur(self):
        return self.a.vecteur(self.b)

On peut s'en servir pour

  1. Avoir une représentation paramétrique de la droite;
  2. 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

[modifier | modifier le wikicode]

Le 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):

Équation cartésienne

[modifier | modifier le wikicode]
    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

[modifier | modifier le wikicode]

L'é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

[modifier | modifier le wikicode]

Le 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 | modifier le wikicode]
    def oalo(self):
        return self.a.y-self.cd()*self.a.x

On 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

[modifier | modifier le wikicode]

Pour 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

[modifier | modifier le wikicode]

Pour savoir si deux droites sont parallèles, on peut

  1. Regarder si leurs vecteurs directeurs sont colinéaires (mieux);
  2. 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é

[modifier | modifier le wikicode]

Pour 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

[modifier | modifier le wikicode]

En 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:

  1. forward ou fd: Pour avancer (l'unité de distance est le pixel)
  2. backward ou bk: Pour reculer
  3. left ou lt: Pour tourner à gauche (l'unité d'angle est le degré)
  4. right ou rt: Pour tourner à droite
  5. goto pour téléporter la tortue (donner l'abscisse puis l'ordonnée)
  6. penup() ou pu() pour que les déplacements de la tortue cessent de laisser des traces à l'écran
  7. pendown() ou pd() pour que les déplacements de la tortue recommencent à laisser des traces
  8. position() renvoie les coordonnées de la tortue
  9. home() renvoie la tortue au centre de l'écran (sa position initiale)
  10. reset() fait pareil mais en effaçant l'écran
  11. circle dessine un cercle (le rayon est en pixels)
  12. 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

[modifier | modifier le wikicode]

Les 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

[modifier | modifier le wikicode]

Pour 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))


Pour 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

[modifier | modifier le wikicode]

Une fois qu'on a vu des nombres négatifs, on peut chercher comment réaliser des opérations dessus:

Pour 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:

Pour 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

[modifier | modifier le wikicode]

De 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)


On a vu dans un chapitre précédent comment la module turtle permet de représenter graphiquement une fonction.

Chute d'une bille sur la planche de Galton

[modifier | modifier le wikicode]

La 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

[modifier | modifier le wikicode]

Pour 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

[modifier | modifier le wikicode]

Une 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):


La 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

[modifier | modifier le wikicode]

Comme 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:

  1. A (comme avance) pour fd(100);
  2. p (comme plus) pour lt(60);
  3. 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.

Bundesarchiv Bild 183-1990-0703-023, Dresden, Rauschgiftspürhund im Einsatz

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

[modifier | modifier le wikicode]

Pour 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

[modifier | modifier le wikicode]

Pour 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):

Sur ce site, il y a un wikibook sur LOGO dont beaucoup d'idées sont transposables ici.


Résolution de systèmes en Python

Un petit exercice
Au village de Trokhafairtrade dans le Swazibwana occidental, on pratique un mélange de commerce équitable et de troc. Ainsi, un habitant a acquis 2 youkoulélés d'Hawaii contre 3 xylophones et 1 € et un écotouriste a acquis un xylophone et un youkoulélé pour la modique somme de 8 €. On demande le prix, en euros, d'un xylophone et le prix d'un youkoulélé sur le marché de Trokhafairtrade.

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

[modifier | modifier le wikicode]

Une 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

[modifier | modifier le wikicode]

Si 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

[modifier | modifier le wikicode]

La 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)

Le 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é).

Les 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 :

Pour 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:

Énoncé

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

[modifier | modifier le wikicode]

Sujet 2005, Exercice 2

[modifier | modifier le wikicode]
Énoncé
Pour organiser une grande manifestation sportive, le professeur d’éducation physique doit rassembler sur le stade un important groupe d’élèves. Le nombre d’élèves est compris entre 2 800 et 2 900. Il en profite pour leur faire remarquer que, regroupés par 2, puis par 3, puis par 4, puis par 5, puis par 6, il en reste toujours 1 ; mais, ô miracle, en se regroupant par 7, il ne reste personne.

On 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

[modifier | modifier le wikicode]

On 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

[modifier | modifier le wikicode]

Maintenant 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

[modifier | modifier le wikicode]

Cette 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!

Maintenant 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!

On 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

[modifier | modifier le wikicode]

Autant 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

[modifier | modifier le wikicode]
Énoncé
Chaque semaine, Jean ramasse entre 40 et 200 œufs qu’il va vendre au marché.

Ce 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.

On 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!

s10=[n for n in s6 if n%10==2]

print(s10)

Plus que 5 nombres à examiner!

solutions=[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.

En voici l'énoncé traduit du Néerlandais à l'Anglais puis au Français :
On choisit au hasard deux entiers x et y strictement supérieurs à 1, x étant le plus petit des deux, et on donne à Sam leur somme qui est inférieure à 100, et à Polly leur produit. Après un temps de réflexion suffisamment long, le dialogue suivant se déroule entre Sam et Polly :
  • 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.

Pour 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

[modifier | modifier le wikicode]

L'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

[modifier | modifier le wikicode]

Si 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

[modifier | modifier le wikicode]

La 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

[modifier | modifier le wikicode]

Puisque 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

[modifier | modifier le wikicode]

Trouver 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 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.