La résolution de systèmes linéaires \(Ax = b\) est l’un des problèmes fondamentaux de l’algèbre linéaire numérique vus en première année de prépa. En première année, on voit généralement que ‘numpy.linalg.solve’ résout ce type de système en une ligne. Mais comprendre ce qui se passe derrière cette fonction opaque est une compétence directement évaluable en concours : comment triangularise-t-on une matrice ? Pourquoi le choix du pivot est-il crucial ? etc. Dans cet article, nous allons implémenter la méthode du pivot de Gauss, étudier l’importance du pivot partiel afin de perfectionner notre première méthode, puis présenter la décomposition LU.
L’élimination de Gauss
Principe
Le principe a déjà été vu en première année de prépa et doit absolument être maîtrisé pour les écrits, quelles que soient les écoles préparées. Résumons rapidement la logique.
L’élimination de Gauss transforme le système \(Ax = b\) en un système triangulaire supérieur équivalent, que l’on résout ensuite par substitution arrière. On travaille sur la matrice augmentée \([A|b]\) de taille \(n \times (n+1)\).
Pour chaque colonne \(k\) de \(0\) à \(n-2\), on utilise la ligne \(k\) comme ligne pivot et on élimine les entrées sous la diagonale : pour chaque ligne \(i > k\), on soustrait \(\frac{a_{ik}}{a_{kk}}\) fois la ligne \(k\) à la ligne \(i\). Le coefficient \(a_{kk}\) est appelé ici pivot, comme tu l’auras sûrement déjà vu.
Implémentation
Voici une manière d’implémenter manuellement cette méthode en Python :
Analysons précisément ce script
On commence par construire la matrice augmentée \(M\) en collant \(b\) comme dernière colonne de \(A\) avec np.hstack. La boucle externe sur \(k\) parcourt chaque colonne pivot de gauche à droite. Pour chaque pivot \(k\), la boucle interne sur \(i\) parcourt toutes les lignes situées en dessous. On calcule le facteur d’élimination \(M[i, k] / M[k, k]\), puis on soustrait ce facteur multiplié par la ligne pivot à la ligne \(i\), ce qui annule le coefficient en position \((i, k)\). On réitère le processus sur chaque pivot avec ‘for k in range’.
Une fois la triangularisation terminée, on initialise le vecteur solution \(x\) à zéro et on remonte de la dernière ligne à la première, exactement de la même manière qu’on le ferait à la main. Pour chaque ligne \(i\), les inconnues \(x[i+1], …, x[n-1]\) sont déjà calculées : \(np.dot(M[i, i+1:n], x[i+1:n])\) calcule leur contribution, on la soustrait du second membre \(M[i, -1]\), et on divise par le pivot diagonal \(M[i, i]\) pour obtenir \(x[i]\). La fonction renvoie le vecteur solution \(x\).
On voit donc que cette fonction permet de résoudre le système linéaire avec la méthode du pivot de Gauss, comme on le ferait manuellement !
Voici un exemple avec ce script :
Avec ces matrices, on obtiendra comme solution : [ 12.5 –11.57142857 5.42857143]
Autrement dit, le système suivant a comme solution \(
\left(\frac{25}{2},-\frac{81}{7},\frac{38}{7}\right)
\) :
\[
\begin{cases}
2x+y-z=8,\\
4x+5y+2z=3,\\
-2x+y+6z=-4.
\end{cases}
\]
Complexité
Lorsqu’on programme une méthode numérique, il ne suffit pas qu’elle fournisse le bon résultat. Effectivement, il convient de savoir combien d’opérations elle effectue lorsque la taille du problème augmente afin d’optimiser l’algorithme pour qu’il soit toujours efficace lorsqu’on manipule des objets de grande dimension.
Dire qu’un algorithme est de complexité (O(n)) signifie que son temps d’exécution croît proportionnellement à (n), tandis qu’une complexité \(O(n^2)\) indique une croissance quadratique, c’est-à-dire que lorsque la taille de la matrice est multipliée par 10, le nombre d’opérations est alors multiplié approximativement par 100.
Ici, la boucle double d’élimination effectue de l’ordre de \(\frac{n^3}{3}\) opérations arithmétiques. La complexité totale est donc en \(O(n^3)\). C’est également la complexité de ‘numpy.linalg.solve’, qui utilise en interne la même approche (ainsi c’est rassurant de trouver la même complexité) !
Le pivot partiel
Le problème du pivot nul ou petit
L’algorithme de Gauss naïf (vu précédemment) échoue si un pivot \(a_{kk}\) est nul, puisque dans le script, la division par zéro est impossible. Mais il peut aussi donner des résultats très imprécis si le pivot est petit en valeur absolue. En effet, le facteur \(\frac{a_{ik}}{a_{kk}}\) peut alors être très grand, ce qui amplifie les erreurs d’arrondi et propage des imprécisions dans tout le système.
Voici un exemple concret : le système \(\begin{pmatrix} 10^{-15} & 1 \\ 1 & 1 \end{pmatrix} x = \begin{pmatrix} 1 \\ 2 \end{pmatrix}\) a pour solution exacte \(x \approx (1, 1)\), mais Gauss sans pivot peut donner une solution très éloignée à cause des erreurs d’arrondi.
La stratégie du pivot partiel
La stratégie du pivot partiel consiste, à chaque étape \(k\), à chercher parmi les lignes \(k, k+1, \dots, n-1\) celle dont le coefficient en colonne \(k\) est le plus grand en valeur absolue et à l’échanger avec la ligne \(k\) avant d’effectuer l’élimination. Cela garantit que le pivot est toujours le plus grand disponible, ce qui minimise l’amplification des erreurs. Présentons donc une version améliorée du pivot de Gauss.
Implémentation avec pivot partiel
L’algorithme reprend le principe du pivot de Gauss classique, mais ajoute une étape de sélection du pivot avant chaque élimination. Pour une colonne \(k\), l’expression M[k:,k] extrait les coefficients situés sous le pivot potentiel, puis np.abs calcule leurs valeurs absolues afin de comparer leur taille. La fonction np.argmax renvoie alors l’indice du plus grand coefficient en valeur absolue, ce qui permet d’identifier le meilleur pivot disponible. Si ce pivot n’est pas déjà sur la ligne \(k\), les lignes sont échangées grâce à l’instruction M[[k,idx\_max]] = M[[idx\_max,k]]. On vérifie ensuite que le pivot n’est pas trop proche de zéro afin d’éviter les divisions dangereuses.
L’élimination se poursuit alors de manière classique en annulant les coefficients situés sous le pivot à l’aide du facteur \(\frac{a_{ik}}{a_{kk}}\). Une fois la matrice transformée en matrice triangulaire supérieure, la substitution arrière permet de calculer successivement toutes les inconnues. Cette stratégie dite du pivot partiel améliore considérablement la stabilité numérique de l’algorithme et constitue la méthode utilisée dans la plupart des bibliothèques de calcul scientifique.
Comparaison numérique
On obtient alors [0.99920072 1. ] avec la première méthode contre [1. 1.] avec la méthode du pivot de Gauss amélioré. Cette deuxième méthode est donc bien plus satisfaisante, comme tu peux le constater, car elle trouve le résultat exact.
La décomposition LU
Principe
Voici un article sur ce type de décomposition matricielle si tu veux approfondir la notion. Ce qui suit constitue simplement une introduction au propos.
La décomposition LU factorise la matrice \(A\) en un produit \(A = LU\), où \(L\) est une matrice triangulaire inférieure à diagonale unité et \(U\) est une matrice triangulaire supérieure. Une fois cette factorisation calculée, résoudre \(Ax = b\) se ramène à résoudre deux systèmes triangulaires :
\[Ly = b \quad \text{puis} \quad Ux = y\]
chacun résolu en \(O(n^2)\) par substitution. L’intérêt est que, si l’on doit résoudre plusieurs systèmes avec la même matrice \(A\), mais des seconds membres \(b\) différents, on calcule \(LU\) une seule fois en \(O(n^3)\) et chaque nouveau système ne coûte que \(O(n^2)\), ce qui est plus optimal que la méthode du pivot de Gauss dans le cas d’un problème avec des matrices de grandes dimensions.
Implémentation avec scipy
En Python, ‘scipy.linalg.lu’ calcule la décomposition LU avec pivot partiel, renvoyant trois matrices \(P\), \(L\), \(U\) telles que \(PA = LU\), où \(P\) est une matrice de permutation.
Note : dans la pratique, on utilise souvent \(PA = LU\) plutôt que \(A = PLU\).
‘solve_triangular’ résout un système triangulaire sans recalculer de décomposition : c’est plus efficace que d’appeler ‘np.linalg.solve’ à nouveau. Le paramètre ‘lower=True’ indique que la matrice est triangulaire inférieure pour la première résolution, et ‘lower=False’ pour la triangulaire supérieure. On vérifie que \(PA = LU\) en utilisant ‘P @ b’ comme second membre du premier système.
Conclusion
En définitive, la résolution de systèmes linéaires repose sur trois algorithmes complémentaires : l’élimination de Gauss, simple à implémenter et à comprendre ; le pivot partiel, indispensable pour compléter la première approche ; et la décomposition LU, plus efficace dans des cas de grandes dimensions. Savoir implémenter Gauss from scratch, expliquer pourquoi le pivot partiel est nécessaire et comprendre le principe de LU sont des compétences pertinentes pour savoir ce qui se cache derrière les fonctions directement proposées par Python pour résoudre des systèmes linéaires.
Il n’y a plus qu’à croiser les doigts pour qu’un tel thème tombe aux concours !
Tu peux retrouver ici le méga-répertoire qui contient toutes les annales de concours et les corrigés. Tu peux également accéder ici à toutes nos autres ressources mathématiques !
