dimanche 11 février 2018
Rudiments d'algèbre linéaire avec Sympy
. Cet article s'adresse à des gens qui ont peu fait d'algèbre linéaire ou beacoup oublié.
Table des matières
On suppose bien entendu que le paquet python-sympy
est installé. Pour les cas particuliers voyez la documentation.
Ensuite pour l'utiliser il faut naturellement importer d'une manière ou d'une autre le module sympy
.
Par exemple on importe tout :
from sympy import *
Ou on importe avec un alias
import sympy as sp
Dans de cas il faut tout préfixer avec sp
.
Le tutoriel est une excellente introduction(1)
J'utiliserai parfois la nouvelle syntaxe de formatage des chaines : Formated string literals. Exemple :
>>> a, b = 0, 1
>>> print(f'a = {a}, b = {b}')
a = 0, b = 1
Le cœur de l'algèbre linéaire dans Sympy
est dans le module matrices
. Aussi nous allons commencer par celles-ci. Dans la suite on suppose que nous importons le module Sympy
en entier.
>>> from Sympy import *
Une matrice $(n, p)$ ou $n×p$ est un tableau de nombres réels (ou de nombres complexes) disposés en lignes et en colonnes. $n$ est le nombre de lignes, et $p$ le nombre de colonnes. Notations :
$$M= \left[\begin{array}{ccc} m_{11} & ... & m_{1p} \\ ... & ... & ... \\ m_{n1} & .. & m_{np} \end{array}\right]$$
On écrit de façon abrégée $M=(m_{ij})$ où $i\in \{1, ..., n\}$ $j\in \{1, ...,p \}$
Les lignes $[m_{i1} ... m_{ip}]$ sont appelées les vecteurs lignes de $M$. Ce sont des matrices $1×p$.
Les colonnes de $M~:~\left[\begin{array}{c}m_{1j}\\ ... \\m_{nj}\end{array}\right]$ sont les vecteurs colonnes de $M$.
Exemple :
>>> M=Matrix([[1, 2, 3, 4], [5, 6, 7, 8]]) # 2 lignes 3 colonnes
>>> M
Matrix([
[1, 2, 3, 4],
[5, 6, 7, 8]])
>>> pprint(M)
⎡1 2 3 4⎤
⎢ ⎥
⎣5 6 7 8⎦
Remarque : Dans sympy pprint
(pretty print) produit un affichage amélioré
M.rows
et M.cols
contiennent les dimensions de la matrice M
.
Sympy
utilise une indexation étendue M[i, j]
analogue à $m_{ij}$, sauf qu'elle commence en zéro.
Avec l'exemple juste au dessus :
>>> print(', '.join([f'M[0,{j}]={M[0,j]}' for j in range(4)]))
M[0,0]=1, M[0,1]=2, M[0,2]=3, M[0,3]=4
>>> print(', '.join([f'M[1,{j}]={M[0,j]}' for j in range(4)]))
M[1,0]=1, M[1,1]=2, M[1,2]=3, M[1,3]=4
Ont peur créer avec Sympy
des varaibles indicées avec la syntaxe : var('x1:4(1:4)')
qui va créer les variables x11, x12, x13, x21, x22, x23, x31, x32, x33
.
Cf. help(Sympy.var)
et help(Sympy.symbols)
.
>>> v = var('x1:4(1:4)')
>>> X = Matrix(3, 3, v); pprint(X)
⎡x₁₁ x₁₂ x₁₃⎤
⎢ ⎥
⎢x₂₁ x₂₂ x₂₃⎥
⎢ ⎥
⎣x₃₁ x₃₂ x₃₃⎦
La transposée d'une matrice $n×p$ $M=(m_{ij})$ est la matrice $p×n$, notée $M^t$, $M^t=(p_{ij})$ telle que $\forall i\in\{1, ..., n\} \forall j\in\{1, ...,p\} p_{ji}=m_{ij}$. Trivialement les lignes deviennet les colonnes et vice versa.
On note aussi la transposée de $M,~{}^tM$.
Avec Sympy
la transposée de M
s'obtient avec M.T
ou M.transpose()
ou transpose(M)
.
>>> M=Matrix(2,3, list(range(6))); pprint(M)
⎡0 1 2⎤
⎢ ⎥
⎣3 4 5⎦
>>> pprint(M.T)
⎡0 3⎤
⎢ ⎥
⎢1 4⎥
⎢ ⎥
⎣2 5⎦
L'instruction Matrix([1, 2, 3])
produit une matrice (un vecteur) colonne [[1], [2],[3]]
. Le constructeur Matrix a converti la liste [1, 2, 3]
en une liste de 3 lignes, réduites chacune à un seul élément (colonne).
>>> C=Matrix([1, 2, -3]); pprint(C)
⎡1 ⎤
⎢ ⎥
⎢2 ⎥
⎢ ⎥
⎣-3⎦
>>> pprint(C.transpose()) # obtention d'un vecteur ligne
[1 2 -3]
On peut forcer le constructeur à produire une matrice ligne avec Matrix(1, 3, [1, 2, 3])
. Dans cette forme Matrix(n, p, [x1, ..., xn])
n et p sont les dimensions de la matrice et la liste est la liste des termes dans l'ordre des lignes.
En accord avec l'indexation étendue, Sympy
autorise des sections étendues de la forme M[r,:]
ou M[:,c]
pour accéder respectivement à la ligne r
ou à la colonne c
. Exemple :
>>> A=Matrix(2,3, list(range(-2, 4))); pprint(A)
⎡-2 -1 0⎤
⎢ ⎥
⎣1 2 3⎦
>>> pprint(A[0,:]) # ligne 0
[-2 -1 0]
>>> pprint(A[:,2]) # colonne 1
⎡0⎤
⎢ ⎥
⎣3⎦
Comme en python l[:]
désigne toute la liste.
La notation M[1:, 1:3]
est aussi possible. Ici on obtient une matrice extraite :
>>> M=Matrix(3, 3, list(range(1, 10))); pprint(M)
⎡1 2 3⎤
⎢ ⎥
⎢4 5 6⎥
⎢ ⎥
⎣7 8 9⎦
>>> pprint(M[1:, 1:3])
⎡5 6⎤
⎢ ⎥
⎣8 9⎦
En math on est souvent amené à écrire une matrice à partir de ses vecteurs colonnes (matrice de changement de base, matrice d'un endomorphisme, ...). La saisie se fait alors comme avec les lignes, puis on transpose la matrice.
>>> var('(u:w)1:4')
(u1, u2, u3, v1, v2, v3, w1, w2, w3)
>>> M=Matrix(3,3, _).T ; M
Matrix([
[u1, u2, u3],
[v1, v2, v3],
[w1, w2, w3]])
J'utilise le script sym.py, en démarrant l'interpréteur ainsi:
python -i sym.py
Il importe Sympy
en entier et comporte deux fonctions matrix
et cmatrix
(inspirées de matrix de numpy) qui permettent une saisie moins verbeuse. En outre ces fonctions convertissent les franctions rationnelles p/q
en Rational(p, q)
avant que python ne s'en soit emparé et ne les ait converties en décimal(2)
Voir la documentation avec help(matrix)
et help(cmatrix)
. Plus accessoirement ce script définit une classe Vect qui mime notre écriture mathématique des vecteurs (lignes).
Si $n = p$ on dit que la matrice est carrée n.
>>> A=Matrix(3,3, list(range(1, 10))); pprint(A)
⎡1 2 3⎤
⎢ ⎥
⎢4 5 6⎥
⎢ ⎥
⎣7 8 9⎦
Si $A=(a_{ij})$ et $B=(b_{ij}$ sont des matrices de mêmes dimensions, la somme est la matrice $C=(c_{ij})=(a_{ij}+b_{ij})$. Exemple :
$$\left(\begin{array} {cc} 1 & 2 \\ 3 & 4 \end{array} \right) + \left(\begin{array}{cc} -1 & 1 \\ 1 & 4 \end{array} \right) = \left(\begin{array}{cc} 0 & 3 \\ 4 & 8 \end{array} \right)$$
>>> A,B=Matrix([[1, 2],[3, 4]]), Matrix([[-1, 1], [1, 4]])
>>> pprint(A+B)
⎡0 3⎤
⎢ ⎥
⎣4 8⎦
Multiplier une matrice par un nombre, c'est multiplier chaque terme par ce nombre $k.(a_{ij})=(ka_{ij})$. Exemple, avec la matrice A ci-dessus :
>>> A = Matrix([[1, 2],[3, 4]])
>>> var('k') # on introduit la variable k
k
>>> B = k*A
>>> pprint(B)
⎡ k 2⋅k⎤
⎢ ⎥
⎣3⋅k 4⋅k⎦
>>> B.subs(k, -2); B
Matrix([
[-2, -4],
[-6, -8]])
Prenons les matrices carrées d'ordre $n$ (par exemple $n=3$). Désignons par O la matrice bourrée de zéros et par -A la matrice faite avec les opposés des termes de A.
L'addition a les mêmes propriétés que l'addition des nombres :
$$ \begin{array}{ll}(A+B)+C=A+(B+C) & A+B=B+A \\ A+O=O+A=A & A+(-A)=(-A)+A=O\end{array}$$
Par ailleurs la multiplication par un nombre a des propriétés commodes :
$$ \begin{array}{ll} 1.A=A &(\lambda+\mu).A=\lambda.A+\mu.A \\ \lambda.(A+B)=\lambda.A+\lambda.B & \lambda.(\mu.A)=(\lambda \mu).A \end{array}$$
On dit que l'ensemble des matrices carrées d'ordre $n$ est un espace vectoriel sur ℝ. (Cf. infra). Ce résultat se généralise à l'ensemble des matrices de dimensions $(n, p)$ entiers.
Avec Sympy
, la matrice d'ordre n
nulle est zeros(n)
:
>>> zeros(2)
Matrix([
[0, 0],
[0, 0]])
>>> var('a:d')
(a, b, c, d)
>>> M=Matrix([[a, b], [c, d]]); pprint(M)
⎡a b⎤
⎢ ⎥
⎣c d⎦
>>> M.det()
a*d - b*c
Le déterminant de la matrice $M=\left(\begin{array}{cc} a & b \\ c & d \end{array}\right)$ est le nombre $det(M)=ad -bc$.
>>> var('m1:4(1:4)')
(m11, m12, m13, m21, m22, m23, m31, m32, m33)
>>> M = Matrix(3, 3, _); pprint(M)
⎡m₁₁ m₁₂ m₁₃⎤
⎢ ⎥
⎢m₂₁ m₂₂ m₂₃⎥
⎢ ⎥
⎣m₃₁ m₃₂ m₃₃⎦
>>> M.det()
m11*m22*m33 - m11*m23*m32 - m12*m21*m33 + m12*m23*m31 + m13*m21*m32 - m13*m22*m31
Le déterminant de la matrice $M=\left(\begin{array}{ccc} m_{11} & m_{12} & m_{13} \\ m_{21} & m_{22} & m_{23} \\ m_{31} & m_{32} & m_{33} \end{array}\right)$ est le nombre :
$$\det(M)=m_{11}m_{22}m_{33}- m_{11}m_{23}m_{32}- m_{12}m_{21}m_{33}+ m_{12}m_{23}m_{31}+ m_{13}m_{21}m_{32}- m_{13}m_{22}m_{31}$$
Prenons deux matrices carrées $A$ et $B$ de rang $n$. On définit leur produit comme suit : le terme $m_{ij}$ est la somme $m_{ij}=\sum_{k=1}^{n}a_{ik}b_{kj}$. C'est une sorte de produit scalaire du vecteur ligne $i$ de $A$de $A$ par le vecteur colonne $j$ de $B$.
>>> a = var('a1:3(1:3)')
>>> b = var('b1:3(1:3)')
>>> A=Matrix(2, 2, a); pprint(A)
⎡a₁₁ a₁₂⎤
⎢ ⎥
⎣a₂₁ a₂₂⎦
>>> B=Matrix(2, 2, b); pprint(B)
⎡b₁₁ b₁₂⎤
⎢ ⎥
⎣b₂₁ b₂₂⎦
>>> P = A*B; pprint(P)
⎡a₁₁⋅b₁₁ + a₁₂⋅b₂₁ a₁₁⋅b₁₂ + a₁₂⋅b₂₂⎤
⎢ ⎥
⎣a₂₁⋅b₁₁ + a₂₂⋅b₂₁ a₂₁⋅b₁₂ + a₂₂⋅b₂₂⎦
Si on désigne par $1_2$ la matrice $\left(\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right)$ on a alors toujours $1_2A=A1_2=A$. Le nom de la matrice unité dans Sympy
est eye
:
>>> eye(2)
Matrix([
[1, 0],
[0, 1]])
On se limitera aux matrices carrées.
Le produit matriciel est associatif : $(A.B).C=A.(B.C)$
En général il ne commute pas : $A.B \ne B.A$
>>> A, B = Matrix([[1, 2], [4, 1]]), Matrix([[1, 1], [1, 0]])
>>> A*B == B*A
False
La matrice diagonale $I_n=\left( \begin{array}{cccc}1 & 0 & ... & 0 \\ 0 & 1 & ... & 0 \\ . & . & ... & . \\ 0 & 0 & ... & 1 \end{array} \right)$ est l'élément neutre de la multiplication des matrices : $I_n.A=A.I_n=A$.
En général les matrices ne sont pas inversibles.
Une matrice carrée est inversible si et seulement son déterminant n'est pas nul.
>>> M=Matrix(3, 3, list(range(9)))
>>> M.det()
0
>>> M.inv()
Traceback (most recent call last):
...
ValueError: Matrix det == 0; not invertible.
Bien que :
>>> B = Matrix([[1, 1], [1, 0]]); pprint(B)
⎡1 1⎤
⎢ ⎥
⎣1 0⎦
>>> pprint(B.inv())
⎡0 1 ⎤
⎢ ⎥
⎣1 -1⎦
Pour davantage d'explications voir Matrices (linear algebra) et Python numbers vs. Sympy Numbers
Un espace vectoriel sur ℝ (ou ℂ), est un ensemble E d'éléments appelés vecteurs, qui se prêtent à deux opérations l'addition et la multiplication par un nombre réel ou complexe. Les propriétés sont données en bas de page (2)
L'espace usuel $ℝ^3$ de la géométrie élémentaire est un espace vectoriel sur ℝ.
Voici quelques calculs simples avec Sympy
>>> u, v, w = Vect('1/2, -1, 0'), Vect('1, 0, 3'), Vect('1, 1, 1')
>>> print(f'u = {u}, v = {v}, w = {w}')
u = (1/2, -1, 0), v = (1, 0, 3), w = (1, 1, 1)
>>> print(f'u + v = {u+v}, u - v = {u-v}, 2*u - v + w = {2*u-v+w}')
u + v = (3/2, -1, 3), u - v = (-1/2, -1, -3), 2*u - v + w = (1, -1, -2)
Remarque : Pour la nouvelle notation de formatage des chaines voir Formated string literals
L'ensemble des polynomes P de la forme $P(X)= a_2X^2 + a_1X + a_0$ est un espace vectoriel sur ℝ (ou ℂ).
Nous verrons plus bas que l'ensemble des matrices carrées de rang $n$ est un espace vectoriel su ℝ.
Un sous espace d'un espace vectoriel E est une partie E' de E stable pour les deux opérations : $\vec{0} \in E' \:et\: \forall \vec{u} \in E', \: \forall \vec{v} \in E',\: \forall \lambda \in ℝ, \; \lambda .\vec{u}+\vec{v} \in E'$
Toute combinaison linéaire de p vecteurs $\vec{v}_1, ...,\vec{v}_p$ de E', appartient à E' :
$\forall a_1, ..., a_p \in ℝ^p \; a_1\vec{v}_1+ ...+a_p\vec{v}_p \in E'$
Soit F un sous espace vectoriel d'un espace vectoriel E, et une famille de vecteurs $f_n=(\vec{v_1}, \vec{v_2}, ..., \vec{v_n})$.
$f_n$ est une famille est génératrice de $F$, signifie que tout vecteur $\vec{v}$ de $F$ est une combinaision linéaire des vecteurs de $f_n$ :
$$\forall \vec{v} \in F,\: \exists (α_1, ..., α_n) \in ℝ^n,\: \vec{v}=α_1\vec{v_1}+ ... +α_n\vec{v_n}$$
$F$ est alors le plus petit espace vectoriel engendré par $f_n$.
$f_n$ est libre (ou que ses vecteurs sont linéairement indépedants) si:
$$\exists (α_1, ..., α_n) \in ℝ^n $non tous nuls$,\: α_1\vec{v_1}+ ... +α_n\vec{v_n}=\vec{0}$$
Dans le cas contraire on dit que $f_n$ est liée.
$f_n$ est une base de $F$ si $f_n$ est libre et génératrice. On dit alors que la dimension
de $F$ est $n$.
Une application d'un espace vectoriel $E$ dans un espace vectoriel $F$ est une application qui respecte la linéarité : $f(\lambda \vec{u}+\mu \vec{v})=\lambda \vec{v}+\mu \vec{v}$
Dans le cas particulier $E=F$ on parle d'endomorphisme dans $E$. o
Dans le cas particulier où $F=ℝ$ (ou $F=ℂ$), on parle de forme linéaire.
Ici nous nous limiterons à ces deux cas.
Soit $f$ une application linéaire de $E$ dans $E$ (endomorphisme de $E$) et $(B)=(\vec{e_{1}}, ..., \vec{e_{n}}$ (donc $\dim(E)=n$
Ecrivons la matrice $M=\left(\begin{array}{ccc} m_{11} & ... & m_{1n} \\ ... & ... & ... \\ m_{n1} & .. & m_{nn} \end{array}\right)$ dans laquelle les colonnes sont les images $f(\vec{e_1}), ..., f(\vec{e_n})$ c'est à dire la colonne $j$ est le vecteur $f(\vec{e_j})=m_{1j}\vec{e_1}+ ... + m_{nj}\vec{e_n}$.
Cette matrice est la matrice associée à l'endomorphisme $f$ dans la base $(B)$.
$E$ est le plan géométrique muni de la base $(B)=(\vec{i}, \vec{j})$
Soit $s_1$ et $s_2$ les symétries d'axes $\vec{i}+\vec{j}$ et $\vec{j}$ :
$$s_1(\vec{i})=\vec{j}=0.\vec{i}+1.\vec{j} \:$ et $\: s_1(\vec{j})=\vec{i}=1.\vec{i}+0.\vec{j}$$
$$s_2(\vec{i})=-\vec{i}=-1.\vec{i}+0.\vec{j} \:$ et $\: s_2(\vec{j})=\vec{j}=0.\vec{i}+1.\vec{j}$$
>>> S1 = Matrix(([0, 1], [1, 0])).T; S1 # .T transpose ici met en colonne
Matrix([
[0, 1],
[1, 0]])
>>> S2 = Matrix(([-1, 0], [0, 1])).T; S2
Matrix([
[-1, 0],
[ 0, 1]])
En composant ces deux symétries on devrait obtenir la rotation $r=s_2\circ s_1$ d'angle $\alpha$ :
$$\alpha=2 imes(\vec{i}+\vec{j},\:\vec{j})=\displaystyle{2 imes\frac{\pi}{4}=\frac{\pi}{2}}$$
>>> R = S2*S1; R
Matrix([
[0, -1],
[1, 0]])
>>> Matrix(([cos(a), sin(a)], [-sin(a),cos(a)])).T.subs(a, pi/2)
Matrix([
[0, -1],
[1, 0]])
On vérifie bien que la matrice de $s_2\circ s_1$ est la même que celle de $rot(\dfrac{\pi}{2})$
Sympy
propose trois rotations de $ℝ^3$ de vecteurs unitaires les vecteurs de la base canonique de $ℝ^3$ : rot_axis1
, rot_axis2
et rot_axi3
. L'argument de la fonction est l'angle de la rotation et le résultat est sa matrice. Exemples de rotations d'angle $\dfrac{\pi}{2}$ :
>>> R1 = rot_axis1(pi/2); R1
Matrix([
[1, 0, 0],
[0, 0, 1],
[0, -1, 0]])
>>> R2 = rot_axis2(pi/2); R2
Matrix([
[0, 0, -1],
[0, 1, 0],
[1, 0, 0]])
>>> R3 = rot_axis3(pi/2); R3
Matrix([
[ 0, 1, 0],
[-1, 0, 0],
[ 0, 0, 1]])
A la différence du plan, le produit de deux rotations n'est pas toujours une rotation.
Plutôt que la formule de Rodrigues, nous allons utiliser ici les matrices de passage.
Soit $(\vec{i},\,\vec{j},\,\vec{k})$ la base canonique de $ℝ^3$
Prenons un vecteur unitaire $\vec{p}(\frac{2}{3},\,\frac{2}{3},\,\frac{1}{3})$ et complétons le par les vecteurs $\vec{p}=\vec{n}\land(\vec{j}+\vec{k})$ et $\vec{q}=\vec{n}\land\vec{p}$ pour obtenir une orthonormée.
#!/usr/bin/env python3
from Sympy import *
# base canonique
i, j, k = (eye(3)[:,c] for c in range(3))
# axe de rotation normalisé
n = Matrix((2/S(3), 2/S(3), 1/S(3)))
# complétion de la base : (n, p, q) orthonormée directe
p = n.cross(j + k)
q = n.cross(p)
# matrice de passage de (i, j, k) à (n, p, q) orthogonale directe
P = (n.row_join(p)).row_join(q)
# Matrice de la rot(n, π/2)
R = P*rot_axis1(pi/2)*P.inv()
La mise en page de la documentation de sympy comporte trois colonnes. Si comme moi cette présentation vous énerve vous avez quatre solutions :
1. Ouvrir la documutation sans javascript, mais alors MathJax ne fonctionnera plus et certaines expressions mathématiques apparaitront en code latex.
2. Télécharger la version pdf.
3. Créer une copie locale «zen» de la doc
Voici une copie d'écran.
Pour cela :
Télécharger une une copie de la doc en html et l'adapter ainsi :
Purger le code html de toutes les lignes commençant avec <link rel="stylesheet" href="http://live.sympy.org/static/
ou donner une adresse factice comm http://zob.sympy.org ...
find -type f -name "*.html" -exec sed -i "s/live\.sypmy\.org/zob.sympy.org/g" {} \;
Bidouiller la feuille de style _static/basic.css
pour virer les colonnes de droite et de gauche.
(2)div.sphinxsidebarwrapper {
display: none;
...
}
div.document,
div.documentwrapper,
div.bodywrapper {
margin: 0 !important;
width: 100%;
}
C'est une question délicate qu'on ne trouve pas dans d'autres logiciels de calcul symbolique. sympy définit bien un type rationnel (Rational(p, q)
), mais il se fait courcircuiter par python lorsque on exécute une fonction, la fonction __new__ de Matrix par exemple, qui convertit immédiatement une écriture fractionnaire en float
.
Exemple 1/2
est converti en 0.5
avant que sympy n'ait pris la main. Les solutions sont alors :
+ entrer la forme Rational(1, 2)
+ entrer S(1)/2
où S
est un raccourci de Sympify
, qui veut dire traduire en sympy. Il suffit que l'un des termes de la fraction soit converti.
+ saisir la donnée sous forme de chaine et faire le traitement dans la fonction.
Rappels des propriétes des espaces vectoriels.
$E$ étant un espace vectoriel
$$ \forall \vec{u} \in E, \forall \vec{v} \in E, \; \vec{u} + \vec{v} \in E,\; (\vec{u} + \vec{v}) + \vec{w} = \vec{u} + (\vec{v} + \vec{w}), \;\vec{u} + \vec{v} = \vec{v} + \vec{u}$$
$$ \exists \vec{0} \in E, \: \forall \vec{u} \in E, \; \vec{0}+\vec{u}=\vec{u} \:\: et \:\: \forall \vec{v} \in E, \; \exists -\vec{v} \in E, \; \vec{v}=-\vec{v}=\vec{0}$$
$$\forall \lambda \in ℝ, \, \forall µ \in ℝ \, \forall \vec{u} \in E \, \forall \vec{v} \in E \, :$$
$$1.\vec{v} = \vec{v}, \; (\lambda + \mu).\vec{u}=\lambda . \vec{u} + \mu . \vec{u}, \; \lambda .(\vec{u}+\vec{v})=\lambda .\vec{u}+\lambda .\vec{v}, \; \lambda .(\mu.\vec{u})=(\lambda \mu). \vec{u}$$
A ce stade il n'y a pas de multiplication entre vecteurs.