Chapitre 37 — faire calculer l'ordinateur sur les objets des Maths expertes : complexes avec cmath, matrices avec numpy, graphes et chaînes de Markov.
Naviguez avec les flèches ← → du clavier, ou via le Sommaire.
complex, module cmathnumpy : produit, puissance, inverse, systèmesL'ordinateur calcule vite, fait voir et permet de conjecturer ; mais il calcule en flottants : on compare toujours « à $10^{-10}$ près », jamais avec ==.
complexLe nombre $i$ s'écrit 1j ; $3+4i$ s'écrit 3 + 4j ou complex(3, 4). Pour un complexe z : z.real, z.imag (des flottants), abs(z) $=\lvert z\rvert$ et z.conjugate() $=\overline{z}$.
>>> z = 3 + 4j >>> abs(z), z * z.conjugate() (5.0, (25+0j)) >>> w = 1 - 2j >>> z * w, z / w ((11-2j), (-1+2j)) >>> 1j ** 2 (-1+0j)
j ** 2 provoque NameError : j seul est un nom de variable. Un calcul avec un complexe reste un complexe, et z < w n'a pas de sens ($\mathbb{C}$ n'est pas ordonné).
cmathcmath.phase(z) est l'argument dans $]-\pi,\pi]$, cmath.polar(z) le couple $(\lvert z\rvert,\arg z)$, cmath.rect(r, t) le complexe $re^{it}$ et cmath.sqrt(z) une racine carrée de $z$. Le module math, lui, ignore les complexes.
>>> import cmath >>> cmath.polar(1 + 1j) (1.4142135623730951, 0.7853981633974483) >>> cmath.sqrt(7 - 24j) (4-3j) >>> cmath.exp(1j * cmath.pi) (-1+1.2246467991473532e-16j) >>> cmath.exp(1j * cmath.pi) == -1 False >>> cmath.isclose(cmath.exp(1j * cmath.pi), -1) True
$e^{i\pi}=-1$, mais $\sin\pi$ est calculé à $10^{-16}$ près : on ne teste jamais l'égalité de deux complexes calculés avec ==, on utilise cmath.isclose.
Les racines $n$-ièmes de l'unité sont les $u_k=e^{i\frac{2k\pi}{n}}$, $k\in\llbracket 0,n-1\rrbracket$ ; celles de $Z=\rho e^{i\alpha}$ sont les $z_0u_k$, avec $z_0=\sqrt[n]{\rho}\,e^{i\frac{\alpha}{n}}$.
def racines_unite(n) : L = [] for k in range(n) : L.append(cmath.exp(2j * k * cmath.pi / n)) return L >>> [arrondi(u) for u in racines_unite(4)] [(1+0j), 1j, (-1+0j), -1j] >>> arrondi(sum(racines_unite(5))) 0j
La fonction arrondi(z) du cours arrondit les deux coordonnées à $10^{-10}$ : sans elle, $u_1=i$ s'affiche avec une partie réelle de l'ordre de $10^{-17}$. On retrouve $\mathbb{U}_4=\{1,i,-1,-i\}$ et la somme nulle des racines.
Si $\delta$ est une racine carrée quelconque de $\Delta=b^2-4ac$ (fournie par cmath.sqrt), les solutions de $az^2+bz+c=0$ sont :
def second_degre(a, b, c) : d = cmath.sqrt(b ** 2 - 4 * a * c) return (-b + d) / (2 * a), (-b - d) / (2 * a) >>> second_degre(1, 2 + 1j, -1 + 7j) ((1-2j), (-3+1j)) >>> second_degre(1, -2, 5) ((1+2j), (1-2j))
Version condensée de la fonction du cours. Le choix de $\delta$ ne change que l'ordre des solutions ; on vérifie $z_1+z_2=-\dfrac{b}{a}$ et $z_1z_2=\dfrac{c}{a}$.
$P$ est codé par la liste de ses cœfficients par degré croissant : P[k] $=a_k$. On calcule de l'intérieur vers l'extérieur :
Soit $n$ multiplications, au lieu de $\dfrac{n(n+1)}{2}$ pour le calcul direct.
def horner(P, z) : n = len(P) - 1 s = P[n] for k in range(n - 1, -1, -1) : s = s * z + P[k] return s >>> horner([-1, 0, 0, 1], 1j) # P(z) = z^3 - 1 (-1-1j) >>> horner([-1 + 7j, 2 + 1j, 1], 1 - 2j) 0j
À partir de $z_0$, on itère $z_{k+1}=z_k-\dfrac{P(z_k)}{P'(z_k)}$ ; rien n'exige que $z_k$ soit réel. $P'$ a pour cœfficients $[a_1,2a_2,\ldots,na_n]$.
def derivee(P) : return [k * P[k] for k in range(1, len(P))] def newton(P, z0, n) : Q = derivee(P) z = z0 for k in range(n) : z = z - horner(P, z) / horner(Q, z) return z >>> newton([-1, 0, 0, 1], 1j, 10) (-0.5+0.8660254037844387j)
La boucle for garantit la terminaison, pas la convergence : on vérifie abs(horner(P, z)) < 1e-10. La racine atteinte (ici $j$) dépend de $z_0$ ; si $P'(z_k)=0$, ZeroDivisionError.
numpynp.array transforme une liste de lignes en matrice. A.shape donne la taille, A[i, j] le cœfficient, A[:, j] une colonne, A.T la transposée ; np.eye(n) $=I_n$.
>>> import numpy as np >>> A = np.array([[2, 1, 0], [-1, 3, 1]]) >>> A.shape (2, 3) >>> A[0, 1] np.int64(1) >>> A.T array([[ 2, -1], [ 1, 3], [ 0, 1]])
Le cœfficient $a_{ij}$ du cours est A[i - 1, j - 1] : source d'erreurs permanente.
@A + B, k * A : somme et produit par un réel. A @ B est le produit matriciel $AB$, défini si le nombre de colonnes de $A$ égale le nombre de lignes de $B$ (sinon ValueError).
>>> A = np.array([[-1, 1], [1, -1]]) >>> B = np.array([[1, -1], [1, -1]]) >>> A @ B array([[0, 0], [0, 0]]) >>> B @ A array([[-2, 2], [-2, 2]]) >>> A * B array([[-1, -1], [ 1, 1]])
A * B multiplie terme à terme et A ** 2 élève chaque cœfficient au carré : ce ne sont ni $AB$ ni $A^2$. On voit aussi $AB\neq BA$ et $AB=0$ avec $A\neq0$, $B\neq0$.
np.linalg.matrix_power(A, n) $=A^n$, np.linalg.det(A) $=\det A$, np.linalg.inv(A) $=A^{-1}$ (sinon LinAlgError) et np.linalg.solve(A, B) résout $AX=B$ pour $A$ inversible.
>>> np.linalg.matrix_power(np.array([[1, 2], [0, 1]]), 10) array([[ 1, 20], [ 0, 1]]) >>> np.linalg.inv(np.array([[1, -1], [1, 3]])) array([[ 0.75, 0.25], [-0.25, 0.25]]) >>> A = np.array([[5, 9], [1, 2]]) >>> np.linalg.solve(A, np.array([3, -2])) array([ 24., -13.])
On retrouve $A^{10}$ avec $A^n=\begin{pmatrix}1&2n\\0&1\end{pmatrix}$, puis l'inverse du cours et la solution $(24,-13)$ du système $\begin{cases}5x+9y=3\\x+2y=-2\end{cases}$.
np.allclose et contre-exemplesAvec $A=\begin{pmatrix}2&1\\1&3\end{pmatrix}$, det renvoie 5.000000000000001 et $AA^{-1}$ contient $-5{,}6\times10^{-17}$ au lieu de $0$ : le test == échoue. On compare deux matrices calculées avec np.allclose.
>>> A = np.array([[1, -1], [1, 3]]) >>> B = np.array([[2, 1], [1, 3]]) >>> inv = np.linalg.inv >>> np.allclose(inv(A @ B), inv(B) @ inv(A)) True >>> C = A @ A + 2 * A @ B + B @ B >>> np.allclose((A + B) @ (A + B), C) False
Un False est un contre-exemple : $(A+B)^2=A^2+2AB+B^2$ est fausse en général. Un True, même répété, ne prouve rien : c'est une conjecture à démontrer.
Chaque bloc de deux lettres (A $=0$, …, Z $=25$) forme une colonne $X$, chiffrée en $Y=KX$ modulo $26$. Pour déchiffrer, $K'\equiv u\begin{pmatrix}d&-b\\-c&a\end{pmatrix}\ [26]$, où $u$ est l'inverse de $\det K$ modulo $26$ : il existe si, et seulement si, $\det K$ et $26$ sont premiers entre eux.
>>> K = np.array([[3, 3], [2, 5]]) >>> pow(9, -1, 26) # det K = 9 3 >>> hill("MATHSX", K) 'KYAVTV' >>> hill("KYAVTV", inverse_mod26(K)) 'MATHSX'
A % 26 réduit chaque cœfficient dans $\llbracket 0,25\rrbracket$, même négatif ; ord("M") - 65 donne $12$ et chr(12 + 65) redonne "M".
Les termes se calculent par une boucle ; l'état stable, solution de $(I-A)S=B$, s'obtient avec np.linalg.solve. Une colonne se code par un tableau à une dimension.
def suite(A, B, X0, n) : X = X0 for k in range(n) : X = A @ X + B return X >>> A = np.array([[0.5, 0.2], [0.1, 0.6]]) >>> B = np.array([1, 2]) >>> suite(A, B, np.array([10, 5]), 50) array([4.44444446, 6.11111113]) >>> np.linalg.solve(np.eye(2) - A, B) array([4.44444444, 6.11111111])
La suite semble converger vers $S=\begin{pmatrix}40/9\\55/9\end{pmatrix}$, cohérent avec $X_n=A^n(X_0-S)+S$ : les cœfficients de $A^{50}$ sont de l'ordre de $10^{-8}$. L'ordinateur suggère, la démonstration reste à faire.
M.sum(axis=1) : sommes des lignes, soit les degrés (sortants si orienté) ; M.sum() // 2 : nombre d'arêtes.np.array_equal(M, M.T) : graphe non orienté ?[i, j] de $M^k$ compte les chaînes de longueur $k$ de i à j.>>> M.sum(axis=1) # graphe G2 du cours array([3, 2, 2, 1, 0]) >>> np.linalg.matrix_power(N, 2) # graphe orienté G11 array([[0, 1, 2, 1, 1], [0, 0, 1, 0, 1], [0, 0, 0, 0, 0], [0, 1, 0, 0, 0], [0, 0, 1, 1, 0]])
Somme des degrés $8=2\times4$ arêtes. Le cœfficient [0, 2] de $N^2$ vaut $2$ : deux chaînes de longueur $2$ du sommet $1$ au sommet $3$ du cours.
On tient une file des sommets découverts : on retire le premier, on ajoute ses voisins non encore vus, et l'on recommence tant que la file n'est pas vide.
def composante(M, s) : n = len(M) vus = [s] file = [s] while file != [] : i = file.pop(0) for j in range(n) : if M[i, j] > 0 and j not in vus : vus.append(j) file.append(j) return vus def connexe(M) : return len(composante(M, 0)) == len(M) >>> composante(M, 0), connexe(M) ([0, 1, 2, 3], False)
Chaque sommet entre au plus une fois dans la file : au plus $n$ tours, donc l'algorithme termine. Le sommet $5$ du cours (4) est isolé : $G_2$ n'est pas connexe.
$\pi_n=\pi_0T^{\,n}$ s'écrit pi0 @ np.linalg.matrix_power(T, n) : un tableau à une dimension placé à gauche de @ est traité comme une matrice ligne.
>>> T = np.array([[0.4, 0.6], [0.8, 0.2]]) >>> pi0 = np.array([1, 0]) >>> pi0 @ T @ T array([0.64, 0.36]) >>> pi0 @ np.linalg.matrix_power(T, 50) array([0.57142857, 0.42857143]) >>> np.array([0, 1]) @ np.linalg.matrix_power(T, 50) array([0.57142857, 0.42857143])
On retrouve $\pi_2$ du cours et la convergence vers $\begin{pmatrix}\frac47&\frac37\end{pmatrix}$, quelle que soit la distribution initiale.
$\pi T=\pi$ équivaut à $(T-I)^{T}\pi^{T}=0$. Une équation de ce système homogène est redondante : on la remplace par $\pi_1+\cdots+\pi_n=1$, puis on résout avec np.linalg.solve.
def invariante(T) : n = len(T) A = (T - np.eye(n)).T A[n - 1] = np.ones(n) # derniere ligne : somme = 1 B = np.zeros(n) B[n - 1] = 1 return np.linalg.solve(A, B) >>> pi = invariante(T) >>> pi, np.allclose(pi @ T, pi) (array([0.57142857, 0.42857143]), True)
Depuis l'état $i$, on tire $r$ dans $[0,1[$ et on choisit l'état $j$ tel que $t_{i0}+\cdots+t_{i,j-1}\leqslant r<t_{i0}+\cdots+t_{ij}$ : il sort avec la probabilité $t_{ij}$.
def etat_suivant(T, i) : r = random.random() j = 0 s = T[i, 0] while r > s and j < len(T) - 1 : j = j + 1 s = s + T[i, j] return j >>> random.seed(2026) >>> simulation(T, 0, 10) [0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0] >>> traj = simulation(T, 0, 10000) >>> traj.count(0) / len(traj) 0.5754424557544245
simulation enchaîne $n$ appels à etat_suivant. La fréquence de l'état 0 approche $\frac47\simeq0{,}571$ : une illustration, pas une démonstration.
1j, z.real, abs(z), z.conjugate()phase, polar, rect, sqrt, iscloseA @ B ; matrix_power(A, n)det, inv, solve(A, B)np.allclose, jamais ==M.sum(axis=1) ; $M^k$ compte les chaînespi0 @ matrix_power(T, n)j seul n'est pas $i$ : on écrit 1j.A * B et A ** 2 calculent terme à terme : le produit matriciel est A @ B.A[i - 1, j - 1].cmath.isclose, np.allclose, jamais == ; un déterminant de $10^{-16}$ est nul.True sur un exemple ne démontre rien ; un False est un contre-exemple.cmath.polar(1 + 1j) ?A ?np.linalg.matrix_power(A, 5), et surtout pas A ** 5.A @ np.linalg.inv(A) vaut $I_2$ ?np.allclose(A @ np.linalg.inv(A), np.eye(2)).[i, j] de $M^3$ ?i au sommet j.pi0 @ np.linalg.matrix_power(T, 10).