6 Algèbre
Tout d’abord une méthode de résolution de système linéaire par un algorithme direct est étudiée, puis une méthode itérative sera utilisée pour déterminer le vecteur propre associé à la plus grande valeur propre d’une matrice. Enfin les groupes générés par un ensemble de permutations seront étudiés.
Exercice 6.1 - Décomposition LU
La décomposition LU consiste à décomposer une matrice \(A\) de taille \(n \times n\) sous la forme \(A=LU\) où \(L\) est une matrice triangulaire inférieure avec des 1 sur la diagonale et \(U\) une matrice triangulaire supérieure. Une fois la décomposition \(A=LU\) d’une matrice connue, il est alors très facile de résoudre le problème linéaire \(A\boldsymbol{x} = \boldsymbol{b}\) en résolvant d’abord \(L\boldsymbol{y} = \boldsymbol{b}\) puis \(U\boldsymbol{x} = \boldsymbol{y}\). L’avantage de la décomposition LU sur l’algorithme de Gauss, par exemple, est qu’une fois la décomposition LU connue, il est possible de résoudre le système linéaire rapidement avec des seconds membres différents.
Vu que \(l_{ik}=0\) si \(k>i\), nous avons:
\[ a_{ij} = \sum_{k=1}^n l_{ik}u_{kj} = l_{ii}u_{ij} + \sum_{k=1}^{i-1} l_{ik}u_{kj} \,, \]
et donc comme \(l_{ii}=1\):
\[ u_{ij} = a_{ij} - \sum_{k=1}^{i-1} l_{ik}u_{kj} \,. \]
D’un autre côté, vu que \(u_{kj}=0\) si \(k>j\), alors:
\[ a_{ij} = \sum_{k=1}^n l_{ik}u_{kj} = l_{ij}u_{jj} + \sum_{k=1}^{j-1} l_{ik}u_{kj} \,, \]
et donc si \(u_{jj}\neq0\):
\[ l_{ij} = \frac{1}{u_{jj}} \left( a_{ij} - \sum_{k=1}^{j-1} l_{ik}u_{kj} \right) \,. \]
Ainsi, si les \((i-1)\) premières lignes de \(U\) et les \((i-1)\) premières colonnes de \(L\) sont connues, il est possible de déterminer la \(i\)-ème ligne de \(U\) par:
\[ u_{ij} = a_{ij} - \sum_{k=1}^{i-1} l_{ik}u_{kj} \,, \quad j \geq i \,, \]
puis la \(i\)-ème colonne de \(L\) par:
\[ l_{ji} = \frac{1}{u_{ii}} \left( a_{ji} - \sum_{k=1}^{i-1} l_{jk}u_{ki} \right) \,, \quad j > i\ \,. \]
Pour que la décomposition LU d’une matrice \(A\) soit possible il faut que les \(u_{ii}\) ne soient jamais nuls. Cela est effectivement le cas lorsque la matrice \(A\) et toutes ses sous-matrices principales sont inversibles.
a. Écrire une fonction LU(A) qui retourne le résultat de la décomposition LU d’une matrice.
b. Donnée la décomposition LU d’une matrice \(A\), écrire une fonction solve(L,U,b) qui résout le système linéaire \(A\boldsymbol{x} = \boldsymbol{b}\).
c. Écrire une fonction LU_inplace(A) qui ne crée pas de nouveaux tableaux pour \(L\) et \(U\) mais modifie \(A\) de sorte que sa partie triangulaire inférieure (sans la diagonale) corresponde à \(L\) et sa partie triangulaire supérieure (avec la diagonale) corresponde à \(U\). Faire également une version solve_inplace prenant en entrée la sortie de LU_inplace et retournant la solution \(\boldsymbol{x}\), sans utiliser de nouveaux tableaux.
d. En utilisant la décomposition LU de la matrice \(A\), écrire une fonction det(A) qui retourne son déterminant.
Exercice 6.2 - Méthode de la puissance itérée
Le but de cet exercice est de déterminer le vecteur propre d’une matrice associé à la valeur propre de plus grand module, en supposant que celle-ci est unique. Étant donné une matrice réelle \(A\) de taille \(n\times n\) et un vecteur \(\boldsymbol{x}_0\in\mathbb{R}^n\), la suite de vecteurs \((\boldsymbol{x}_k)_ {k\in\mathbb{N}}\) est définie par:
\[ \boldsymbol{x}_{k+1} = \frac{A\boldsymbol{x}_k}{\Vert A\boldsymbol{x}_k\Vert} \,. \]
En supposant par exemple que la matrice \(A\) soit diagonalisable, alors il est possible de montrer que la suite \((\boldsymbol{x}_k)_ {k\in\mathbb{N}}\) converge à un signe près vers le vecteur propre associé à la plus grande valeur propre de \(A\) en valeur absolue. La convergence a lieu presque sûrement pour tous les choix de \(\boldsymbol{x}_0\).
a. Définir une fonction puissance(A, x0, k) qui retourne \(\boldsymbol{x}_k\). Avec cette fonction, déterminer le plus grand vecteur propre de la matrice:
\[A = \begin{pmatrix}0.5 & 0.5\\ 0.2 & 0.8 \end{pmatrix} \,.\]
b. Déterminer la valeur propre associée au vecteur propre précédent.
c. Écrire un programme permettant de calculer automatiquement la valeur propre de plus grand module et le vecteur propre associé d’une matrice carrée avec une certaine précision donnée. On choisira le vecteur initial \(\boldsymbol{x}_0\) aléatoirement.
d. En supposant la matrice \(A\) diagonalisable avec une unique valeur propre de plus grand module, montrer que la suite \(\boldsymbol{x}_k\) converge à un signe près vers le vecteur propre associé à cette valeur propre de plus grand module pour presque tous les choix de \(\boldsymbol{x}_0\).
e. Regarder la documentation de NumPy pour trouver les fonctions permettant de calculer les vecteurs propres et valeurs propres d’une matrice.
f. Comparer la rapidité du code précédent et des fonctions NumPy pour des matrices de tailles \(n\times n\) avec \(n=10,100,1 000\).
Exercice 6.3 - Exponentielle de matrices
Le but de cet exercice est de développer un algorithme permettant de calculer l’exponentielle d’une matrice réelle carrée. Si \(A\) est une matrice réelle carrée, son exponentielle est définie par la série:
\[ \exp(A) = \sum_{k=0}^\infty \frac{A^k}{k!} \,, \]
par analogie avec l’exponentielle sur les nombres réels. Ici \(A^k\) représente le produit matriciel de \(A\) avec elle-même \(k\) fois.
a. Définir les tableaux NumPy correspondant aux matrices \(A_1\) et \(A_2\) définies par:
\[ \begin{align*} A_{1} &= \begin{pmatrix} 1 & 0.8 & 0.6\\ 0.8 & 0.2 & 0.8\\ 0 & 1.2 & 0.9 \end{pmatrix} , & A_{2} &= \begin{pmatrix} 2 & 3 & 2\\ 1 & 2 & 3\\ 4 & 3 & 5.2 \end{pmatrix} . \end{align*} \]
b. Définir une fonction matrix_power(A,n=20) retournant une approximation de \(\exp(A)\) obtenue en gardant uniquement les \(n+1\) premiers termes de la série, i.e. les termes de \(k=0\) à \(k=n\).
c. Tester sur les matrices \(A_1\) et \(A_2\) et comparer avec les résultats de la fonction expm du module linalg de Scipy.
d. Sans utiliser la fonction norm de NumPy ou Scipy, définir une fonction calculant la norme de Frobenius \(\Vert A\Vert_{F}\) d’une matrice \(A\) de taille \(m \times m\) définie par:
\[ \Vert A\Vert_{F}^{2}=\operatorname{tr}(AA^{t})=\sum_{i=1}^{m}\sum_{j=1}^{m}|a_{ij}|^{2} \,. \]
Calculer les normes de Frobenius des matrices \(A_1\) et \(A_2\).
e. Pour les matrices \(A_1\) et \(A_2\), tracer l’erreur en norme de Frobenius entre le résultat de matrix_power et le résultat de expm en fonction du nombre de termes \(n\) gardés. Mettre une échelle logarithmique sur l’axe des ordonnées.
D’un point de vue théorique, il est possible de montrer que l’erreur est bornée par:
\[ \left\Vert \exp(A)-\sum_{k=0}^{n}\frac{A^{k}}{k!}\right\Vert_{F} \leq\frac{e^{\Vert A\Vert_{F}}}{(n+1)!}\Vert A\Vert_{F}^{n+1} \,. \]
f. Tracer cette borne en fonction de \(n\) pour différentes valeurs de la norme \(\Vert A\Vert_{F}\) allant de \(2\) à \(20\), avec également une échelle logarithmique sur les ordonnées. En déduire grossièrement le nombre de termes à garder pour que la borne soit inférieure à la précision machine de \(10^{-15}\) si \(\Vert A\Vert_{F}=20\). Comparer la borne théorique avec ce qui a été observé à la question précédente.
Une idée basique pour améliorer la convergence de la série lorsque la norme de la matrice est grande est d’effectuer un changement d’échelle à l’aide de la relation:
\[ \exp{A} = \bigl(\exp(A/p)\bigr)^p \,. \]
pour \(p \geq 1\) un entier bien choisi tel que \(\Vert A/p \Vert_F\) soit petite, par exemple inférieure à un.
g. En utilisant la propriété précédente, écrire une fonction matrix_power_opt(A,n=20) fondée sur cette propriété.
h. Refaire le même graphique qu’au point e mais avec cette nouvelle fonction et commenter.
Exercice 6.4 - Groupes de permutations
Le but est d’étudier les groupes de permutations en développant des algorithmes pour les caractériser. Un groupe de permutations \(G \subset S_n\) est défini comme étant généré par un certain nombre de permutations: \(G = \langle g_1, g_2,\dots,g_k \rangle\), avec \(g_i\in S_n\) une permutation de l’ensemble \(\{1,2,\dots,n\}\). Une permutation
\[g = \begin{pmatrix}1 & 2 & 3 & 4\\ 4 & 1 & 2 & 3 \end{pmatrix} ,\]
peut être représentée en Python par le tuple g = (0, 4, 1, 2, 3). Le zéro est ajouté afin de se conformer avec l’indexation à partir de zéro de Python et ainsi de simplifier un peu les implémentations. Cela veut simplement dire que le sommet 0 va sur le sommet 0. À noter que cet exercice se prête particulièrement bien à une implémentation orientée objet, et dans ce cas les questions peuvent être adaptées en conséquence.
a. Écrire une fonction product(g1,g2) qui calcule le produit de deux permutations.
b. Écrire une fonction inverse(g) qui calcule l’inverse d’une permutation.
c. Écrire une fonction qui retourne la décomposition d’une permutation sous forme de cycles représentés par une liste du tuples.
d. Écrire une fonction qui retourne la permutation correspondant à une liste de cycles, c’est-à-dire qui fait l’inverse de la fonction précédente.
e. ! En Python un groupe \(G = \langle g_1,g_2,\dots,g_k \rangle\) engendré par une famille de permutations, peut être représenté par une liste de permutations G = [g1,g2,...,gk]. Écrire une fonction orbit(G,x) qui retourne l’orbite d’un point \(x\in\{1,2,\dots,n\}\):
\[ O_x = Gx = \{gx, g \in G\} \,. \]
f. !! Définir une fonction stabilizer(G,x) qui retourne le stabilisateur d’un point \(x\in\{1,2,\dots,n\}\):
\[ G_x = \{g \in G : g x = x \} \,, \]
sous la forme d’un ensemble de générateurs.