Modélisation mathématique#
Sage est orienté objet#
Python et Sage utilisent fortement la programmation orientée objet.
Même si cela reste relativement transparent pour une utilisation
occasionnelle, il est utile d’en savoir un minimum, d’autant que ce fait
est très naturel dans un contexte mathématique.
Le paradigme de la programmation orientée objet s’appuie sur un
principe: «toute entité du monde physique ou mathématique que l’on
souhaite manipuler avec l’ordinateur est modélisé par un objet»; de
plus cet objet est une instance d’une classe. Ainsi, le nombre
rationnel \(o=12/35\) est modélisé par un objet qui est une instance de la
classe Rational :
12/35
o = 12/35
type(o)
Noter que cette classe est vraiment associée à l’objet \(12/35\), et non
seulement à la variable o qui le contient:
type(12/35)
Précisons les définitions. Un objet est une portion de la mémoire de l’ordinateur qui contient l’information nécessaire pour représenter l’entité qu’il modélise. La classe quant à elle définit deux choses:
la structure de données d’un objet, c’est-à-dire comment l’information est organisée dans le bloc mémoire. Par exemple, la classe
Rationalstipule qu’un nombre rationel comme \(12/35\) est représenté, en gros, par deux nombres entiers: son numérateur et son dénominateur.son comportement, et en particulier les opérations sur cet objet: par exemple comment on extrait le numérateur d’un nombre rationel, comment on calcule sa valeur absolue, comment on multiplie ou additionne deux nombres rationels, etc. Chacune de ces opération est implantée par une méthode (respectivement
numer,abs, {__mult__}, {__add__}, …).
Pour factoriser un nombre entier \(o\), on va donc appeller la méthode
factor avec la syntaxe suivante:
o = 720
o.factor()
que l’on peut lire comme: «prendre la valeur de o et lui appliquer la
méthode factor sans autre argument». Sous le capot, effectue le calcul
suivant:
type(o).factor(o)
De gauche à droite: «demander à la classe de (la valeur de) o
(type(o)) la méthode appropriée de factorisation (type(o).factor),
et l’appliquer à o».
Notons au passage que l’on peut appliquer une opération à une valeur sans passer par une variable:
720.factor()
et donc en particulier enchaîner les opérations, de la gauche vers la droite. Ici, on prend le numérateur d’un nombre rationnel, que l’on factorise ensuite:
o = 720 / 133
o.numerator().factor()
Applications#
Polymorphisme#
En quoi cela nous concerne-t-il? Tout d’abord, l’orientation objet
permet le polymorphisme: quelque soit l’objet o que l’on veut
factoriser, on peut toujours utiliser la notation o.factor() (ou son
raccourci factor(o)). De même, calquant les notations mathématiques
usuelles, le produit de deux objets a et b peut toujours être noté
a*b même si l’algorithme utilisé dans chaque cas est différent (Pour
une opération arithmétique binaire comme le produit, la procédure de
sélection de la méthode appropriée est un peu plus compliquée que ce qui
a été décrit précédemment. En effet elle doit gérer des opérations
mixtes comme la somme \(2 + 3/4\) d’un entier et d’un nombre rationnel. En
l’occurence, \(2\) sera converti en nombre rationnel \(2/1\) avant
l’addition. C’est le modèle de coercion de Sage qui est en charge de
cela.). Voici un produit de deux nombres entiers:
3 * 7
un produit de deux nombres rationnels, obtenu par produit des numérateurs et dénominateurs puis réduction:
(2/3) * (6/5)
Un produit de deux nombres complexes, utilisant \(I^2=-1\) :
(1 + I) * (1 - I)
des produits commutatifs formels de deux expressions:
(x + 2) * (x + 1)
(x + 1) * (x + 2)
Outre la simplicité de notation, cela permet d’écrire des programmes génériques comme:
def puissance_quatre(a):
a = a * a
a = a * a
return a
qui s’appliquent à tout objet admettant les opérations utilisées (ici la multiplication):
puissance_quatre(2)
puissance_quatre(3/2)
puissance_quatre(I)
puissance_quatre(x+1)
M = matrix([[0,-1],[1,0]]); M
puissance_quatre(M)
Introspection#
Plus prosaïquement, l’orientation objet permet l’introspection: on peut ainsi accéder à l’aide en ligne spécifique à la factorisation des nombres entiers avec:
o = 720
o.factor?
voire à l’implantation de cette méthode, précédée de son aide en ligne:
o.factor??
En passant au dessus des détails techniques, on distingue bien que
Sage délègue le calcul à d’autres logiciels (Pari, Kash, …).
Enfin, on peut utiliser la complétion automatique pour demander
interactivement à un objet o quelles sont toutes les méthodes que
l’on peut lui appliquer:
o.<tab>
Ici, il y en a beaucoup; voici celles dont le nom commencent par n :
o.n<tab>
Voici toutes les méthodes dont le nom contient prime :
o.*prime*?
Éléments, parents, catégories#
Éléments et parents#
Dans la section précédente, nous avons vu la notion technique de
classe d’un objet. Dans la pratique, il est suffisant de savoir que
cette notion existe; on a rarement besoin de regarder explicitement le
type d’un objet. En revanche Sage introduit une contrepartie plus
conceptuelle de cette notion que nous allons aborder maintenant: celle
du parent d’un objet.
Supposons par exemple que l’on veuille déterminer si un élément \(a\) est inversible. La réponse ne va pas seulement dépendre de l’élément lui-même, mais de l’ensemble \(A\) auquel il est considéré appartenir. Par exemple, le nombre \(5\) n’est pas inversible dans l’ensemble \(\ZZ\) des entiers, son inverse \(1/5\) n’étant pas un entier:
a = 5; a
a.is_unit()
En revanche, il est inversible dans l’ensemble des rationnels:
a = 5/1; a
a.is_unit()
Comme nous l’avons vu dans la section précédente, Sage répond
différemment à ces deux questions car les éléments \(5\) et \(5/1\) sont
dans des classes différentes:
type(5)
type(5/1)
Dans certains systèmes de calcul formel orientés objet, tels que MuPAD
ou Axiom l’ensemble \(X\) auquel \(x\) est considéré appartenir (ici \(\ZZ\)
ou \(\QQ\)) est simplement modélisé par la classe de \(x\). Sage suit
l’approche de Magma, et modélise \(X\) par un objet supplémentaire
associé à \(x\), et appelé son parent:
parent(5)
parent(5/1)
On peut retrouver ces deux ensembles avec les raccourcis:
ZZ
QQ
et les utiliser pour convertir aisément un élément de l’un à l’autre lorsque cela a un sens:
QQ(5).parent()
ZZ(5/1).parent()
ZZ(1/5)
Voici \(1\) en tant qu’entier \(1\in\ZZ\), nombre rationnel \(1\in\QQ\), et approximations flottantes \(1{,}0\in\RR\) et \(1{,}0+0{,}0i \in\CC\) :
ZZ(1), QQ(1), RR(1), CC(1)
Exemple: Combinatoire#
Selon le même principe, lorsque l’on demande toutes les partitions de l’entier 5, le résultat est un objet qui modélise cet ensemble:
P = Partitions(5); P
Pour obtenir la liste de ces objets, il faut le demander explicitement:
P.list()
Cela permet de manipuler formellement des grands ensembles:
Partitions(100000).cardinality()
Et de calculer paresseusement avec. Ici, on tire au hasard une main de cinq cartes à jouer:
Symboles = Set(["Coeur", "Carreau", "Pique", "Trefle"])
Valeurs = Set([2, 3, 4, 5, 6, 7, 8, 9, 10, "Valet", "Dame", "Roi", "As"])
Cartes = cartesian_product([Valeurs, Symboles])
Mains = Subsets(Cartes, 5)
Mains.cardinality()
Mains.random_element() # random
et là on manipule un mot infini défini comme point fixe d’un morphisme:
m = WordMorphism('a->acabb,b->bcacacbb,c->baba')
m.fixed_point('a')
Complément: Constructions#
Les parents étant eux-même des objets, on peut leur appliquer des opérations. Ainsi, on peut construire le produit cartésien \(\QQ^2\) :
cartesian_product([QQ, QQ])
retrouver \(\QQ\) comme corps des fractions de \(\ZZ\) :
ZZ.fraction_field()
construire l’anneau des polynômes en \(x\) à coefficients dans \(\ZZ\) :
ZZ['x']
Par empilement successif, on peut construire des structure algébriques avancées comme l’espace des matrices \(3\times 3\) à coefficients polynomiaux sur un corps fini:
Z5 = GF(5); Z5
P = Z5['x']; P
M = MatrixSpace(P, 3, 3); M
dont voici un élément:
m = M.random_element(); # random
m
m.det()
Exemple: algèbre linéaire#
Dans les exemples ci-dessous, nous ferons de l’algèbre linéaire sur le corps fini \(\ZZ/7\ZZ\) :
K = GF(7); K
list(K)
Nous avons choisi ce corps à titre d’illustration pour avoir des résultats lisibles. On aurait pu prendre des coefficients entiers, rationnels, ou numériques à plus ou moins haute précision. Les aspects numériques seront abordés plus en détail dans l’exposé suivant. Notons au passage que, même en calcul exact, il est possible de manipuler de relativement grosses matrices:
n = 500
M = random_matrix(K, n, sparse=True, density=3/n)
M.visualize_structure() # not tested
n = 10000
M = random_matrix(K, n, sparse=True, density=3/n)
M.rank() # random
Définissons donc une matrice à coefficients dans \(\ZZ/7\ZZ\) :
A = matrix(K, 4, [5,5,4,3,0,3,3,4,0,1,5,4,6,0,6,3]); A
Calculons le polynôme caractéristique de cette matrice:
P = A.characteristic_polynomial(); P
On vérifie le théorème de Cayley-Hamilton sur cet exemple:
P(A)
Notons que l’information sur le corps de base est préservée:
P.parent()
ce qui influe directement sur la factorisation de ce polynôme:
factor(P)
factor(x^4 + 5*x^3 + 6*x + 2)
Le calcul ci-dessus nous donne les valeurs propres: -3=4,-6=1 et -5=2. Quels sont les espaces propres?
A.eigenspaces_left()
Récupérons ces espaces propres:
E = dict(A.eigenspaces_left())
E[2]
E[2] n’est pas une liste de vecteurs ni une matrice, mais un objet
qui modélise l’espace propre \(E_2\), comme le sous-espace de
\((\ZZ/7\ZZ)^4\) décrit par sa base échelon réduite. On peut donc lui
poser des questions:
E[2].dimension()
E[2].basis()
V = E[2].ambient_vector_space(); V
Voire faire des calculs avec:
E[2] + E[4]
v = V([1,2,0,3])
v in E[2]
E[2].echelon_coordinates(v)
E[2].is_subspace(E[4])
E[2].is_subspace(V)
Q = V/E[2]; Q
Q( V([0,0,0,1]) )
On veut maintenant manipuler \(A\) comme un morphisme sur \(V\) :
phi = End(V)(A); phi
v = V.an_element()
v
phi(v)
(phi^-1)(v)
phi^4 + 5*phi^3 + 6*phi + 2
(phi - 1).image()
(phi - 1).kernel() == E[1]
phi.restrict(E[2])
En résumé#
« Mathematics is the art of reducing any problem to linear algebra » William Stein
Il serait en principe suffisant d’implanter l’algèbre linéaire sur les matrices
Le pari de Sage: modéliser au plus près les mathématiques, pour que l’utilisateur ou le programmeur puisse s’exprimer dans le langage adapté au problème considéré.
Complément: Catégories#
Un parent n’a, en général, pas lui-même un parent, mais une catégorie qui indique ses propriétés:
C = QQ.category(); C
De fait Sage sait que \(\QQ\) est un corps:
QQ in Fields()
et donc, par exemple, un groupe additif commutatif:
QQ in CommutativeAdditiveGroups()
Voici tous les axiomes satisfaits par \(\QQ\) :
C.axioms()
et les catégories de \(\QQ\) :
G = C.category_graph()
G.set_latex_options(format="dot2tex")
view(G, tightpage=True, viewer="pdf")
Sage en déduit que \(\QQ[x]\) est un anneau euclidien:
QQ['x'].category()
En général, il peut combiner des axiomes et des structures:
Magmas().Associative() & Magmas().Unital().Inverse() & Sets().Finite()
Et appliquer par exemple le théorème de Wedderburn:
Rings().Division() & Sets().Finite()
Toutes ces propriétés sont utilisées pour calculer rigoureusement et plus efficacement sur les éléments de ces ensembles.
Expressions versus domaines de calcul explicites#
Dans cette section, nous donnons quelques exemples typiques pour lesquels il est important de contrôler le domaine de calcul. En première lecture, on peut passer rapidement sur les exemples plus avancés pour arriver directement à la synthèse de fin de section.
Exemple: simplification d’expressions#
Soit \(c\) une expression un tout petit peu compliquée:
a = var('a')
c = (a+1)^2 - (a^2+2*a+1)
et cherchons à résoudre l’équation en \(x\) donnée par \(cx=0\) :
eq = c * x == 0
L’utilisateur imprudent pourrait être tenté de simplifier cette équation par \(c\) avant de la résoudre:
eq2 = eq / c; eq2
solve(eq2, x)
Heureusement, Sage ne fait pas cette erreur:
solve(eq, x)
Ici, Sage a pu résoudre correctement le système car le coefficient \(c\)
est une expression polynomiale. Il est donc facile de tester si \(c\) est
nul; il suffit de le développer:
expand(c)
Et d’utiliser le fait que deux polynômes sous forme développée identiques sont égaux. On dit que la forme développée d’un polynôme est une forme normale.
En revanche, sur un exemple à peine plus compliqué, Sage commet une
erreur:
c = cos(a)^2 + sin(a)^2 - 1
eq = c*x == 0
solve(eq, x)
alors même qu’il sait faire la simplification et même le test à zéro correctement:
c.simplify_trig()
c.is_zero()
Cet exemple illustre l’importance du test de nullité, et plus généralement des formes normales, dans un domaine de calcul. Sans lui, tout calcul faisant intervenir une division devient hasardeux. Les algorithmes comme le pivot de Gauss en algèbre linéaire sont particulièrement sensibles à ces considérations.
Exemples: polynômes et formes normales#
Construisons l’anneau \(\QQ[x_1,x_2,x_3,x_4]\) des polynômes en \(4\) variables:
R = QQ['x1,x2,x3,x4']; R
x1, x2, x3, x4 = R.gens()
Les éléments de \(R\) sont automatiquement représentés sous forme développée:
x1 * (x2 - x3)
qui comme nous l’avons vu est une forme normale. On dit alors que \(R\) est à représentation normale. En particulier le test à zéro y est immédiat:
(x1+x2)*(x1-x2) - (x1^2 -x2^2)
Mais ce n’est pas toujours un avantage. Par exemple, si l’on construit le déterminant de Vandermonde \(\prod_{1\leq i < j \leq n} (x_i-x_j)\) :
prod( (a-b) for (a,b) in Subsets([x1,x2,x3,x4],2) )
on obtient \(4!=24\) termes. Alors que la même construction avec une expression reste sous forme factorisée qui est ici beaucoup plus compacte et lisible:
x1, x2, x3, x4 = var('x1, x2, x3, x4')
prod( (a-b) for (a,b) in Subsets([x1,x2,x3,x4],2) )
De même, une représentation factorisée ou partiellement factorisée permet des calculs de { pgcd} bien plus rapides. Réciproquement, il ne serait pas judicieux de mettre automatiquement tout polynôme sous forme factorisée, même s’il s’agit aussi d’une forme normale, car la factorisation est coûteuse et non compatible avec l’addition.
De manière générale, selon le type de calcul voulu, la représentation idéale d’un élément n’est pas toujours sa forme normale. Cela amène les systèmes de calcul formel à un compromis avec les expressions. Un certain nombre de simplifications basiques, comme la réduction des rationnels ou la multiplication par zéro, y sont effectuées automatiquement; les autres transformations sont laissées à l’initiative de l’utilisateur auquel des commandes spécialisées sont proposées.
Exemple: factorisation des polynômes#
Considérons la factorisation de l’expression polynomiale suivante:
x = var('x')
p = 54*x^4+36*x^3-102*x^2-72*x-12
factor(p)
Cette réponse est-elle satisfaisante? Il s’agit bien d’une factorisation
de \(p\), mais son optimalité dépend fortement du contexte! Pour le moment
Sage considère p comme une expression symbolique, qui se trouve être
polynomiale. Il ne peut pas savoir si l’on souhaite factoriser \(p\) en
tant que produit de polynômes à coefficients entiers ou à coefficients
rationnels (par exemple). Pour prendre le contrôle, nous allons préciser
dans quel ensemble (domaine de calcul?) nous souhaitons considérer \(p\).
Pour commencer, nous allons considérer \(p\) comme un polynôme à
coefficient entiers. Nous définissons donc l’anneau \(R=\ZZ[x]\) de ces
polynômes:
R = ZZ['x']; R
Puis nous convertissons \(p\) dans cet anneau:
q = R(p); q
À l’affichage on ne voit pas de différence, mais \(q\) sait qu’il est un élément de \(R\) :
parent(q)
Du coup, sa factorisation est sans ambiguïté:
factor(q)
On procède de même sur le corps des rationels:
R = QQ['x']; R
q = R(p); q
factor(R(p))
Dans ce nouveau contexte, la factorisation est encore non ambiguë; mais
différente de précédemment. Notons au passage que Sage sait que \(R\)
est un anneau euclidien:
R.category()
et donc en particulier un anneau où la factorisation est unique (voir Figure {fig:premierspas:catégories}).
Cherchons maintenant une factorisation complète sur les nombres complexes. Une première option est de s’autoriser une approximation numérique des nombres complexes avec 16 bits de précision:
R = ComplexField()['x']; R
q = R(p); q
factor(R(p))
Une autre est d’agrandir un peu le corps des rationnels; ici, on va rajouter \(\sqrt{2}\).
R = QQ[sqrt(2)]['x']; R
q = R(p); q
factor(R(p))
Enfin, peut-être souhaite-t’on que les coefficients soient considérés modulo \(5\)?
R = GF(5)['x']; R
q = R(p); q
factor(R(p))
Synthèse#
Dans les exemples précédents, nous avons illustré comment l’utilisateur peut contrôler le niveau de rigueur dans ses calculs. D’un côté il peut utiliser les expressions symboliques. Ces expressions vivent dans l’anneau des expressions symboliques:
parent(sin(x))
que l’on peut aussi obtenir avec:
SR
Les propriétés de cet anneau sont assez floues; il est commutatif:
SR.category()
et les règles de calcul font en gros l’hypothèse que toutes les
variables symboliques sont à valeur dans \(\CC\). Le domaine de calcul
(expressions polynomiale? rationnelles? trigonométriques?) n’étant pas
spécifié explicitement, le résultat d’un calcul nécessite le plus
souvent des transformations manuelles pour être mis sous la forme
désirée (voir {sec:calculus:simplifications}), en utilisant par exemple
expand, combine, collect et simplify. Pour bien utiliser ces
fonctions, il faut savoir quel type de transformations elles effectuent
et à quel domaine de calcul ces transformations s’appliquent. Ainsi,
l’usage aveugle de la fonction simplify peut conduire à des résultats
faux. Des variantes de simplify permettent alors de préciser la
simplification à effectuer.
D’un autre côté, l’utilisateur peut construire un parent qui va spécifier explicitement le domaine de calcul. Cela est particulièrement intéressant lorsque ce parent est à forme normale: c’est-à-dire que deux objets éléments sont mathématiquement égaux si et seulement si ils ont la même représentation.
Pour résumer, la souplesse est l’avantage principal des expressions:
pas de déclaration explicite du domaine de calcul;
ajout au vol de nouvelles variables ou fonctions symboliques;
changement au vol du domaine de calcul (par exemple lorsque l’on prend le sinus d’une expression polynomiale);
utilisation de toute la gamme des outils d’analyse (intégration, etc.).
Les avantages de la déclaration explicite du domaine de calcul sont:
vertus pégagogiques: réfléchir au préalable à l’univers où vivent les objets;
rigueur: les résultats obtenus sont garantis corrects (
Sagen’est pas un système de calcul certifié; il peut donc toujours y avoir un bogue informatique; mais il n’y aura pas d’utilisation d’hypothèse implicite);mise sous forme normale automatique (le plus souvent) — cela peut aussi être un inconvénient ! — ;
constructions avancées qui seraient délicates avec des expressions (calculs sur un corps fini ou une extension algébrique de \(\QQ\), dans un anneau non commutatif…).