TP 5 : Modules mathématiques, partie 2 : Networkx et Sympy

1 - Les graphes avec Networkx

Introduction

Mathématiquement, un graphe est donné par un ensemble $S$ dont les éléments s'appellent les sommets, et par un ensemble $A$ dont les éléments sont des paires $\{u, v\}$, avec $u, v\in S$ et $u\ne v$, que l'on appelle les arêtes. (Pour nous, $S$ sera en général fini et composé de nombres entiers.)

On peut représenter un graphe en faisant un dessin simple : un point pour chaque sommet, et un segment pour chaque arête. Le dessin nous aide fortement à réfléchir aux propriétés du graphe.

Il y a plusieurs modules en Python pour travailler avec les graphes, et nous allons voir networkx:

In [1]:
import networkx as nx

Essayons un graphe tout bête :

In [2]:
G= nx.Graph()  # on forme un graphe, qui est vide au début
In [3]:
G.add_nodes_from([0,1,2,3])  # on a ajouté les sommets 0, 1, 2, 3
# notons que pour networkx la traduction de "sommet" est "node", 
# mais c'est souvent aussi "vertex" (pluriel "vertices")
In [4]:
# plaçons quelques arêtes (=edge en anglais):
G.add_edge(0,1)
G.add_edge(1,2)
G.add_edge(2,0)
G.add_edge(0,3)

Le module networkx est capable de faire un petit dessin, mais...

In [5]:
nx.draw(G)

... ça ne marche pas du premier coup! en effet, il faut importer matplotlib. (On se demande pourquoi networkx ne l'importe pas tout seul, mais c'est comme ça.)

In [6]:
import matplotlib
In [7]:
nx.draw(G, with_labels= True)

Le dessin comporte une part de hasard. Essayez with_labels= False pour voir.

Le but des graphes n'est certainement pas de produire des images. Ce qui va nous intéresser, c'est qu'il y a vraiment beaucoup d'algorithmes sur les graphes qui existent déjà, et qui ont été implémentés par exemple dans networkx. Nous allons pouvoir les utiliser à des fins très diverses, même pour des problèmes qui ne semblent pas immédiatement faire intervenir des graphes.

Voyons un exemple. Apparemment, il s'agissait à l'origine d'un défi lancé par un journal à ses lecteurs, mais j'ai perdu la référence :

Quelle est la taille maximale d'un ensemble d'entiers, pris entre 1 et 70, tel que pour toute paire $(a,b)$ de cet ensemble (avec $a\ne b$), la différence $a - b$ n'est pas un carré?

On n'a pas tout de suite l'impression d'un problème sur les graphes. Mais si on pose $S = \{ 1, ..., 70\}$, et si on place une arête entre $i$ et $j$, deux éléments de $S$ avec $i < j$, si $j-i$ n'est pas un carré, alors nous avons un graphe $G$.

Ce que nous cherchons alors, c'est un ensemble de sommets qui soient tous reliés 2 à 2, de taille maximale. (Prenez le temps de bien digérer ça.)

Or (et c'est toujours comme ça!), le concept existe déjà, et les algorithmes aussi : un tel ensemble est appelé une clique du graphe (un mot français qui existe aussi en anglais). Et networkx peut nous trouver les cliques maximales (c'est-à-dire les cliques qu'on ne peut pas agrandir, mais c'est malheureusement différent des cliques de taille maximale, il va nous rester un tout petit peu de travail).

In [8]:
N= 70 # on pourra tester différentes variantes
G= nx.Graph()
G.add_nodes_from(range(1, N+1))

Les carrés qui peuvent intervenir dans les calculs, puisqu'ils sont de la forme $j-i$ avec $i$ et $j$ inférieurs à $N$, doivent eux mêmes être $\le \sqrt N$ :

In [9]:
from math import sqrt as racine
carres= [x**2 for x in range(int(racine(N))) ]

Plaçons alors les arêtes.

In [10]:
for i in range(1, N+1):
    for j in range(i+1, N+1):
        if j - i not in carres:
            G.add_edge(i,j)

On peut essayer de dessiner $G$ mais le résultat est décevant :

In [11]:
nx.draw(G, with_labels= True)

On va maintenant passer en revue les cliques maximales -- données par le module networkx, qui fait le gros du travail pour nous. On en garde une de taille maximale, et c'est fini.

In [12]:
m= 0           # m est la taille de la plus grande clique que nous ayons trouvée jusqu'à présent
solution= []   # solution est la plus grande clique que nous ayons trouvée jusqu'à présent

for clique in nx.find_cliques(G):  # find_cliques donne la liste des cliques maximales
    if len(clique) > m:
        m= len(clique)
        solution= clique
In [13]:
solution
Out[13]:
[1, 6, 45, 11, 21, 8, 13, 58, 63, 66, 53, 16, 51, 40, 3, 68, 48, 34]
In [14]:
pl= nx.draw(G.subgraph(solution))

Exercice ($N=11$). Pour $N=11$, reprendre le code ci-dessus, afficher toutes les cliques maximales, et vérifier qu'elles n'ont pas toutes le même nombre d'éléments.

Quelques commentaires.

A vrai dire, le problème originel était posé pour $N=100$, et pas $N=70$. Avec le code ci-dessus, on n'obtient pas une réponse dans un temps raisonnable. Il existe une bibliothèque de fonctions écrites en C, appelée Cliquer, qui donne une solution presque instantanément :

https://users.aalto.fi/~pat/cliquer.html

C'est peut-être l'occasion de faire de la publicité pour Sagemath. Il s'agit d'une distribution Python avec une grosse collection de bibliothèques/modules à vocation mathématique, dont plusieurs ne sont pas écrits en Python, par exemple Cliquer. (Alors que Anaconda, c'est juste Python.) Avec Sagemath, on peut écrire en gros le même code que ci-dessus (et en fait, c'est même un peu plus court), et il va faire appel à Cliquer sans qu'on le voie, pour donner une réponse immédiate.

Nous n'avons pas fait cette série de TPs avec Sagemath parce que, d'une part, il est plus compliqué à installer, et d'autre part, parce que nous voulons en rester aux outils que vous êtes sûrs à 100% d'utiliser si vous continuez vos études (matplotlib, numpy, networkx et sympy ci-dessous sont vraiment très utilisés, de nos jours).

Allez voir https://www.sagemath.org/.

Retour sur le puzzle du TP 1

Comme promis dans le TP 1, nous allons reprendre le puzzle et le résoudre sans "bidouiller", en laissant le programme faire toutes les recherches pour nous. On va traduire les choses avec des graphes, et vous allez finir avec networkx.

Commençons par un copier-coller de l'énoncé :

Ce puzzle est paru en 1994 dans un magazine de supplément offert avec le Weekend Australian.

Dan est en charge de la promenade de 5 chiens, qui lui ont été confiés par 5 femmes. Seulement, il ne se rappelle plus quel chien appartient à qui, et de plus, il mélange les prénoms des propriétaires. Voici ce dont il est certain.

Le chien de Mme Brown, Loopsie ou Mooksie, n'est pas un terrier. Woopsie, le caniche, n'appartient pas à Mary. Le setter de Marion n'est pas Loopsie. Le basset de Margie s'appelle Smooksie ou Mooksie. Hilary, qui n'est pas Mme Grey, est la maîtresse de Poopsie ou Smooksie. Martha Black n'est pas la maîtresse du dalmatien. Poopsie n'est pas à Mme Grey, mais Mooksie est à Mme White, dont le prénom n'est pas Marion.

Question : sachant que Mme Green est l'une des propriétaires de chien, trouver son prénom, le nom de son chien, et la race de son chien.

(Cette fois-ci, plus d'indication.)

Nous commençons comme dans le TP 1.

In [15]:
from itertools import product

prenoms= ["Mary", "Marion", "Margie", "Hilary", "Martha"]
noms= ["Brown", "Grey", "Black", "White", "Green"]
races= ["terrier", "caniche", "dalmatien", "setter", "basset"]
chiens= ["Loopsie", "Mooksie", "Smooksie", "Poopsie", "Woopsie"]

comb= [ [p, n, r, c] for p,n,r,c in product(prenoms, noms, races, chiens)]

Comme nous avons appris un peu plus de syntaxe entretemps, on peut être un peu plus élégant:

In [16]:
implique=  lambda A, B : B or not A
equivalent=lambda A, B : A == B

On va utiliser aussi la commande all, qui prend une liste de booléens (True/False) et nous dit s'ils sont tous à True :

In [17]:
all([True, False, True])
Out[17]:
False
In [18]:
all([True, True])
Out[18]:
True
In [19]:
filtre= [ [p,n,r,c] for p,n,r,c in comb if all([
    implique(n=="Brown", c in ["Loopsie", "Mooksie"]),
    implique(n=="Brown", r!= "terrier"),
    equivalent(c == "Woopsie", r == "caniche"),
    implique(c=="Woopsie", p != "Mary"),
    equivalent(r=="setter", p=="Marion"),
    implique(p=="Marion", c!="Loopsie"),
    equivalent(r=="basset", p=="Margie"),
    implique(r=="basset", c in ["Smooksie", "Mooksie"]),
    implique(n=="Grey", not(p=="Hilary")),
    implique(p=="Hilary", c in ["Poopsie", "Smooksie"]),
    equivalent(p=="Martha", n=="Black"),
    implique(p=="Martha", not(r=="dalmatien")),
    implique(n=="Grey", not(c=="Poopsie")),
    equivalent(c=="Mooksie", n=="White"),
    implique(n=="White", not(p=="Marion")) ]) # fin du all([...])
    ] # fin du [ de départ
In [20]:
len(filtre)
Out[20]:
27

Nous en étions là, avant de commencer les bidouilles. Rappelons que nous appelons "combinaison" un élément tel que, par exemple:

In [21]:
comb[100]
Out[21]:
['Mary', 'Green', 'terrier', 'Loopsie']

Une solution du puzzle est alors un ensemble de 5 combinaisons (il y a 5 propriétaires de chiens, donc pour chacune il faut lui attribuer un prénom, un chien et une race) qui ne viole aucune des "règles". Or, en plus des règles que nous avons déjà vues, il y en a une qui est implicite : dans une solution, deux combinaisons n'ont aucun caractère en commun (il y a 5 prénoms différents, 5 noms de famille différents, 5 races etc, et il faut les attribuer tous!). Par exemple

['Mary', 'Brown', 'terrier', 'Loopsie']

et

['Mary', 'Green', 'terrier', 'Mooksie']

ne peuvent pas apparaître dans la même solution, puisqu'il n'y a qu'une seule femme prénommée Mary, et par ailleurs un seul des chiens est un terrier.

Mézalors, construisons un graphe $G$ sur les sommets $0, ..., 26$, et plaçons une arête entre $i$ et $j$ si filtre[i] et filtre[j] n'ont aucun caractère en commun. Alors une solution doit former une clique du graphe!

Exercice (puzzle, suite et fin). Finissez le puzzle, en construisant le graphe $G$ et en examinant ses cliques maximales. Vous allez voir au passage que la solution est unique.

2 - Le module généraliste Sympy

Sympy doit son nom au fait qu'il est capable de faire du calcul symbolique, ce qu'on appelle parfois aussi du calcul formel (l'UE que nous faisons s'appelait d'ailleurs il y a quelques années "calcul formel"). On va voir ce que c'est dans un instant, mais il faut tout de suite ajouter que Sympy sait faire bien d'autres choses, et c'est en réalité un module mathématique assez généraliste (même s'il n'empiète pas sur les plates bandes de numpy, matplotlib ou networkx).

Prise en main

Commençons bien sûr par :

In [22]:
import sympy as sp

Et sans attendre, demandons à Sympy de passer par Latex pour afficher tous les résultats (le luxe!) :

In [23]:
sp.init_printing()
In [24]:
sp.sqrt(2)
Out[24]:
$$\sqrt{2}$$

On voit déjà que Sympy possède une fonction sqrt qui est bien différente de celle du module math. Mais on voit aussi que le résultat est joliment formaté (ci-dessous, avec les intégrales et les matrices, ça sera vraiment agréable). C'est la dernière expression d'une cellule qui est traitée de la sorte (si vous utilisez print, tout se passe comme avant).

Si vous voulez afficher plusieurs choses avec ce formatage depuis la même cellule (c'est rare...), on peut le faire avec display. Noter que ce n'est pas sp.display. Par exemple :

In [25]:
for i in range(10):
    display(sp.sqrt(i))
$$0$$
$$1$$
$$\sqrt{2}$$
$$\sqrt{3}$$
$$2$$
$$\sqrt{5}$$
$$\sqrt{6}$$
$$\sqrt{7}$$
$$2 \sqrt{2}$$
$$3$$

Pour commencer, on va constater que Sympy sait faire des choses basiques comme traiter les fractions ou les entiers modulo $n$.

On a déjà vu qu'il y avait plusieurs "modèles" pour les nombres en Python : int et float, mais aussi np.int32 et compagnie... Sympy en introduit plusieurs nouveaux, et celui qui nous intéresse au début est sp.Integer. On s'en sert comme ceci :

In [26]:
x= sp.Integer(3)
x
Out[26]:
$$3$$
In [27]:
type(x)
Out[27]:
sympy.core.numbers.Integer

L'intérêt est qu'un objet créé de la sorte possède beaucoup de méthodes utiles, et par exemple l'opération / marche enfin aussi bien que sur une calculatrice (il était temps) :

In [28]:
x/4
Out[28]:
$$\frac{3}{4}$$

Pour calculer, disons, 7/4 + 3/8, il faut certes taper :

In [29]:
sp.Integer(7)/sp.Integer(4) + sp.Integer(3)/sp.Integer(8)
Out[29]:
$$\frac{17}{8}$$

Bon, nous l'avons vu, Sympy comprend tout seul qu'un objet de type np.Integer divisé par un entier normal doit donner une fraction, donc on peut faire plus simplement :

In [30]:
sp.Integer(7)/4 + sp.Integer(3)/8
Out[30]:
$$\frac{17}{8}$$

Ça reste un peu long. Je propose deux raccourcis ; le premier utilise Z qui évoque $\mathbb{Z}$ :

In [31]:
Z= sp.Integer
In [32]:
Z(7)/4 + Z(3)/8
Out[32]:
$$\frac{17}{8}$$

Autre solution, utiliser:

In [33]:
l= sp.Integer(1)

Il s'agit d'un L minuscule, que l'on ne doit normalement jamais utiliser comme nom de variable, parce que ça ressemble trop à un 1 (l'entier entre 0 et 2). Mais ici, l vaut effectivement 1, et ensuite on peut faire:

In [34]:
7*l/4 + 3*l/8
Out[34]:
$$\frac{17}{8}$$

En effet, 7*l vaut la même chose que sp.Integer(7).

Exercice (le piège des fractions). Expliquez le résultat de sp.Integer(7)/4 + 3/8.

Les objets de type sp.Integer savent faire plein de choses, comme par exemple calculer leur pgcd avec un autre entier :

In [35]:
Z(15).gcd(24)   # gcd= greatest common divisor. Ici on calcule pgcd(15, 24).
Out[35]:
$$3$$

En utilisant Sympy, on va sans le savoir manipuler des objets ayant des types très variés. Par exemple, si on définit :

In [36]:
x= Z(2)

et qu'on prend sa racine carrée :

In [37]:
y= sp.sqrt(x)
y
Out[37]:
$$\sqrt{2}$$
In [38]:
type(y)
Out[38]:
sympy.core.power.Pow

On voit que le type de y est surprenant! Quoique, bien sûr, $\sqrt{2}$ n'est pas un entier, donc la réponse ne pouvait pas être trop simple non plus...

Dans un premier temps, il ne faut pas s'en inquiéter. La seule précaution à prendre, nous l'avons vu ci-dessus, c'est de faire attention quand on mélange les objets Sympy et les objets "normaux" comme les entiers python par exemple. A partir du moment où on construit tous nos objets avec Sympy, on peut lui faire confiance et ne pas chercher à savoir quel est le type de chaque objet.

(Par contre, quand on cherche à devenir vraiment bon avec Sympy, ce n'est plus vrai, évidemment...)

Pour travailler avec des nombres complexes, on peut d'ailleurs faire :

In [39]:
i= sp.sqrt(-1)
In [40]:
i
Out[40]:
$$i$$
In [41]:
i**2
Out[41]:
$$-1$$

Voici des exemples de calculs :

In [42]:
z= 3 + 2*i
z
Out[42]:
$$3 + 2 i$$
In [43]:
z**3
Out[43]:
$$\left(3 + 2 i\right)^{3}$$
In [44]:
(z**3).expand()
Out[44]:
$$-9 + 46 i$$

Nous en disons un peu plus sur la méthode expand ci-dessous.

Avant ça, terminons cette "prise en main" par l'utilisation de Sympy pour travailler avec $\mathbb{Z}/p\mathbb{Z}$ pour $p$ premier -- comme vous le savez sans doute, c'est pour $p$ premier, et seulement dans ce cas, que $\mathbb{Z}/p\mathbb{Z}$ est un corps.

In [45]:
Zmod3= sp.FiniteField(3)  # field en anglais= corps (enfin, en mathématiques)
In [46]:
x= Zmod3(18)
print(x)
0 mod 3

On voit ici que, ayant construit l'objet Zmod3, qui est censé être l'anneau $\mathbb{Z}/3\mathbb{Z}$ lui-même, on peut faire Zmod3(n)n est un entier pour construire un nouvel objet, qui représente n mod 3 bien sûr, et qui possède des méthodes un peu comme celles que nous avons faites nous-mêmes dans le TP 4. Ci-dessus, on a vu que la fonction print fonctionne comme on s'y attend avec un tel entier modulo 3, par exemple. Essayons encore :

In [47]:
p= 11
A= sp.FiniteField(p)  # A comme anneau
for a in range(p):  # voyons encore le petit théorème de Fermat
    print(a, "^", p-1, "=", A(a)**(p-1))
0 ^ 10 = 0 mod 11
1 ^ 10 = 1 mod 11
2 ^ 10 = 1 mod 11
3 ^ 10 = 1 mod 11
4 ^ 10 = 1 mod 11
5 ^ 10 = 1 mod 11
6 ^ 10 = 1 mod 11
7 ^ 10 = 1 mod 11
8 ^ 10 = 1 mod 11
9 ^ 10 = 1 mod 11
10 ^ 10 = 1 mod 11

Ici le but est simplement de vous montrer que les éléments de sp.FiniteField(p) comprennent l'opération puissance **.

(On pourrait lister bien des choses que Sympy sait déjà faire : les nombres binomiaux, le crible d'Eratosthène...)

Le calcul formel

Pour certains, c'est le grand intérêt de l'utilisation des ordinateurs en mathématiques. En réalité, il s'agit d'un ensemble de méthodes pour faire les calculs qu'on vous demande en analyse, en gros.

Commençons par demander à Sympy de nous donner un "symbole" -- on dit parfois une "variable formelle".

In [48]:
x= sp.Symbol("x")  # essayez de changer l'un des deux "x" en "y", et voyez le résultat!
In [49]:
x
Out[49]:
$$x$$

On peut alors construire des expressions plus compliquées (ce sont tous des objets (au sens Python) de type "expression formelle", mais on ne va pas détailler ça).

In [50]:
expr= (x+1)**12
In [51]:
expr
Out[51]:
$$\left(x + 1\right)^{12}$$
In [52]:
expr.expand()
Out[52]:
$$x^{12} + 12 x^{11} + 66 x^{10} + 220 x^{9} + 495 x^{8} + 792 x^{7} + 924 x^{6} + 792 x^{5} + 495 x^{4} + 220 x^{3} + 66 x^{2} + 12 x + 1$$
In [53]:
expr
Out[53]:
$$\left(x + 1\right)^{12}$$

La méthode expand permet de développer -- ou plutôt, de calculer une version développée, mais on voit que l'objet expr n'est pas modifié. On peut aussi faire:

In [54]:
trig= sp.cos(x)**2 + sp.sin(x)**2
In [55]:
trig
Out[55]:
$$\sin^{2}{\left (x \right )} + \cos^{2}{\left (x \right )}$$

Note : bien sûr il a fallu prendre le cosinus et le sinus "version sympy", de même que dans le TP précédent on avait le sinus "version numpy" pour les tableaux.

In [56]:
trig.simplify()
Out[56]:
$$1$$

Et bien sûr la méthode simplify permet de renvoyer une version simplifiée, au gré de ce que Sympy sait faire. Il existe aussi trigsimp pour indiquer à Sympy qu'il faut essayer de la trigo. Dans notre cas les deux marchent:

In [57]:
trig.trigsimp()
Out[57]:
$$1$$

Une chose bien agréable, c'est que Sympy sait dériver et intégrer :

In [58]:
sp.diff(sp.sin(x), x)
Out[58]:
$$\cos{\left (x \right )}$$
In [59]:
sp.diff(sp.exp(x**2), x)
Out[59]:
$$2 x e^{x^{2}}$$
In [60]:
sp.integrate(1/x, x)
Out[60]:
$$\log{\left (x \right )}$$
In [61]:
sp.integrate(x**x, x)
Out[61]:
$$\int x^{x}\, dx$$

Cette dernière réponse, tout à fait correcte, signifie que Sympy ne voit pas comment simplifier l'expression.

Vous ne savez plus par coeur le développement limité de $tan(x)$ en 0 ? Honte à vous, mais vous pouvez le demander à Sympy :

In [62]:
expr= sp.tan(x)
In [63]:
expr.series(x, 0, 12)
Out[63]:
$$x + \frac{x^{3}}{3} + \frac{2 x^{5}}{15} + \frac{17 x^{7}}{315} + \frac{62 x^{9}}{2835} + \frac{1382 x^{11}}{155925} + \mathcal{O}\left(x^{12}\right)$$

Pour récupérer juste la partie polynomiale, c'est un peu compliqué mais c'est possible:

In [64]:
DL= expr.series(x, 0, n= None)  # n= None veut dire que l'ordre n'est pas encore précisé
termes= [ next(DL) for i in range(6)]  # les 6 premiers termes non nuls
P= sum(termes)
P
Out[64]:
$$\frac{1382 x^{11}}{155925} + \frac{62 x^{9}}{2835} + \frac{17 x^{7}}{315} + \frac{2 x^{5}}{15} + \frac{x^{3}}{3} + x$$

(Nous n'expliquerons pas la signification de la commande next, c'est une longue histoire.)

A la rigueur, même si vous ne connaissez plus la formule de Taylor-Young, on peut aussi la demander à Sympy, comme ceci.

In [65]:
f= sp.Function("f")  # ceci dit que f est un "symbole de fonction"
x0= sp.Symbol("x0") # x0 est un symbole comme x
In [66]:
f(x).series(x, x0, 4)  # formule de Taylor-Young pour f en x0 à l'ordre 4
Out[66]:
$$f{\left (x_{0} \right )} + \left(x - x_{0}\right) \left. \frac{d}{d \xi_{1}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \frac{1}{2} \left(x - x_{0}\right)^{2} \left. \frac{d^{2}}{d \xi_{1}^{2}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \frac{1}{6} \left(x - x_{0}\right)^{3} \left. \frac{d^{3}}{d \xi_{1}^{3}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=x_{0} }} + \mathcal{O}\left(\left(x - x_{0}\right)^{4}; x\rightarrow x_{0}\right)$$

Lorsqu'on a une expression formelle, on peut évidemment avoir envie de l'évaluer pour avoir une réponse numérique. On fait comme sur l'exemple suivant.

In [67]:
expr= sp.exp(x**2)+ sp.cos(x)
In [68]:
expr
Out[68]:
$$e^{x^{2}} + \cos{\left (x \right )}$$
In [69]:
expr.subs(x, 0)   # ceci signifie "remplacer x par 0 partout dans l'expression"
Out[69]:
$$2$$

Exercice (graphe d'une expression). Ecrire une fonction trace(expr, x, a, b)expr est une expression qui dépend du symbole x et a et b sont des nombres réels, qui trace le graphe de l'expression sur l'intervalle $[a,b]$ à l'aide de matplotlib (en prenant disons 100 valeurs, cf le TP précédent).

Exercice (développements limités). Ecrire une fonction approx_taylor(expr, x, x0, n, a, b) qui prend une expression expr, un symbole x, un nombre x0, un entier n, et deux nombres a et b, et qui va faire la chose suivante :

  • calculer les $n$ premiers termes non-nuls du DL de expr (qui est une expression qui dépend de x) en x0
  • calculer pour chaque $i \le n$ le polynôme $P_i$ obtenu en ne gardant que les $i$ premiers termes non-nuls
  • tracer sur un même dessin les graphes des $P_i$ pour tout $i$ ainsi que le graphe de expr, sur l'intervalle $[a,b]$, à l'aide de matplotlib, et en s'inspirant d'un exemple du TP 4.

Matrices

Sympy sait faire beaucoup de choses avec les matrices. On peut en construire une comme ceci :

In [70]:
sp.Matrix([ [1,2], [3,4]])
Out[70]:
$$\left[\begin{matrix}1 & 2\\3 & 4\end{matrix}\right]$$

Pour faire la matrice $A=(a_{ij})$ de taille $8\times 8$ avec $a_{ij}= max(i,j)$ (par exemple), on peut faire :

In [71]:
sp.Matrix(8,8, lambda i, j : max(i,j))
Out[71]:
$$\left[\begin{matrix}0 & 1 & 2 & 3 & 4 & 5 & 6 & 7\\1 & 1 & 2 & 3 & 4 & 5 & 6 & 7\\2 & 2 & 2 & 3 & 4 & 5 & 6 & 7\\3 & 3 & 3 & 3 & 4 & 5 & 6 & 7\\4 & 4 & 4 & 4 & 4 & 5 & 6 & 7\\5 & 5 & 5 & 5 & 5 & 5 & 6 & 7\\6 & 6 & 6 & 6 & 6 & 6 & 6 & 7\\7 & 7 & 7 & 7 & 7 & 7 & 7 & 7\end{matrix}\right]$$

La fonction sp.Matrix est donc très flexible, elle peut prendre des arguments très différents -- une liste de listes dans le premier cas, et 2 entiers suivis d'une fonction lambda dans le deuxième cas! (Rares sont les langages de programmation qui permettent ça, et c'est vrai que ça crée souvent la confusion, pour être honnête.)

Pour former rapidement quelques matrices particulières :

In [72]:
# l'identité
sp.eye(4)   # blagounette des auteurs de sympy, car I se prononce comme eye en anglais
Out[72]:
$$\left[\begin{matrix}1 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 1\end{matrix}\right]$$
In [73]:
# la matrice nulle
sp.zeros(4,3)
Out[73]:
$$\left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]$$
In [74]:
# matrice diagonale
sp.diag(1,2,3)
Out[74]:
$$\left[\begin{matrix}1 & 0 & 0\\0 & 2 & 0\\0 & 0 & 3\end{matrix}\right]$$
In [75]:
# matrice diagonale par block
block= sp.Matrix([[1,2], [3,4]])
sp.diag(block, -block, 3*block)
Out[75]:
$$\left[\begin{matrix}1 & 2 & 0 & 0 & 0 & 0\\3 & 4 & 0 & 0 & 0 & 0\\0 & 0 & -1 & -2 & 0 & 0\\0 & 0 & -3 & -4 & 0 & 0\\0 & 0 & 0 & 0 & 3 & 6\\0 & 0 & 0 & 0 & 9 & 12\end{matrix}\right]$$

Tout ce qu'on vous demande de faire avec les matrices en cours (pour ce qui est de la partie calculatoire, en tout cas!) est faisable avec sympy :

In [76]:
A= sp.randMatrix(4)
A
Out[76]:
$$\left[\begin{matrix}5 & 23 & 83 & 13\\71 & 93 & 2 & 15\\84 & 2 & 68 & 11\\73 & 86 & 75 & 84\end{matrix}\right]$$
In [77]:
# determinant
A.det()
Out[77]:
$$-43649969$$
In [78]:
# forme échelonnée
A.echelon_form()
Out[78]:
$$\left[\begin{matrix}5 & 23 & 83 & 13\\0 & -1168 & -5883 & -848\\0 & 0 & -3560950 & -418640\\0 & 0 & 0 & 1274579094800\end{matrix}\right]$$
In [79]:
# inverse (quand il existe)
A**(-1)
Out[79]:
$$\left[\begin{matrix}- \frac{370577}{43649969} & \frac{102901}{43649969} & \frac{474898}{43649969} & - \frac{23213}{43649969}\\\frac{354511}{43649969} & \frac{477196}{43649969} & - \frac{341585}{43649969} & - \frac{95347}{43649969}\\\frac{530597}{43649969} & - \frac{55695}{43649969} & \frac{88256}{43649969} & - \frac{83728}{43649969}\\- \frac{514650}{43649969} & - \frac{528256}{43649969} & - \frac{141791}{43649969} & \frac{712190}{43649969}\end{matrix}\right]$$
In [80]:
# une base du noyau (noyau= ensemble des v tels que A*v = 0)
A= sp.Matrix([ [0,1], [0,1]  ])  # faites-le à la main pour comparer!
A.nullspace()
Out[80]:
$$\left [ \left[\begin{matrix}1\\0\end{matrix}\right]\right ]$$
In [81]:
# resoudre le système linéaire A*v = B où B est une colonne
# prenons A et B:
A= sp.randMatrix(2)
B= sp.Matrix([[1], [2]])
A, B
Out[81]:
$$\left ( \left[\begin{matrix}32 & 62\\54 & 75\end{matrix}\right], \quad \left[\begin{matrix}1\\2\end{matrix}\right]\right )$$
In [82]:
A.solve(B)
Out[82]:
$$\left[\begin{matrix}\frac{49}{948}\\- \frac{5}{474}\end{matrix}\right]$$
In [83]:
# trouver une base de l'espace vectoriel engendré par les colonnes:
A= sp.Matrix([ [1, 2, 3],  [0, 4, 4], [1, 0, 1]  ]) # la 3e colonne est la somme des 2 premières
A.columnspace()
Out[83]:
$$\left [ \left[\begin{matrix}1\\0\\1\end{matrix}\right], \quad \left[\begin{matrix}2\\4\\0\end{matrix}\right]\right ]$$

Exercice (matrice). Prenez une feuille d'exercices sur les matrices de ce semestre, ou du précédent, et faites-en 4 ou 5 à l'aide de Sympy au lieu de les faire à la main.