Résolution d'un système d'équations linéaires

Introduction

Un système linéaire d'équations est un ensemble de deux ou plusieurs équations linéaires impliquant le même ensemble de variables. La solution d'un système linéaire est l'ensemble des valeurs pour les variables qui satisfont simultanément toutes les équations du système.

Par exemple, ceci est un système linéaire d'équations avec deux variables, \(x\) et \(y\): \(\begin{cases} 2x + 3y = 5 \\ 4x - y = 1 \end{cases}\)

Un système linéaire général d'équations peut être écrit comme suit :

Ce système peut être vu comme une équation matricielle \(\mathbf{A} \mathbf{x} = \mathbf{b}\), où :

Sous cette forme, le système linéaire peut être interprété comme une combinaison linéaire de vecteurs. Chaque ligne de la matrice \(\mathbf{A}\) représente les coefficients d'une équation linéaire, et le vecteur \(\mathbf{x}\) contient les variables. Le vecteur du côté droit \(\mathbf{b}\) représente les constantes. Résoudre le système implique de trouver le vecteur \(\mathbf{x}\) tel que lorsque la matrice \(\mathbf{A}\) est multipliée par \(\mathbf{x}\), le résultat est le vecteur \(\mathbf{b}\).

Existence et unicité de la solution

Pour les systèmes linéaires, l'existence et l'unicité des solutions sont bien définies et peuvent être analysées à l'aide de techniques d'algèbre linéaire.

• Le système \( Ax = b \) a une solution si et seulement si la matrice augmentée \([A | b]\) a le même rang que la matrice des coefficients \(A\).
• Si \(\text{rang}(A) = \text{rang}([A | b])\), le système est cohérent (a au moins une solution).
• Si \(\text{rang}(A) = \text{rang}([A | b]) < n\) (où \(n\) est le nombre de variables), le système a une infinité de solutions.
• Si \(\text{rang}(A) < \text{rang}([A | b])\), le système est incohérent (n'a pas de solutions).
• Pour une matrice carrée \(A\), si \(\det(A) \ne 0\) (c'est-à-dire \(A\) est inversible), le système a une solution unique.
• Si \(\det(A) = 0\), le système peut ne pas avoir de solutions ou avoir une infinité de solutions.

Inversion matricielle

Une solution naïve à un système linéaire d'équations \(\mathbf{A} \mathbf{x} = \mathbf{b}\) consiste à trouver l'inverse de la matrice \(\mathbf{A}\) (notée \(\mathbf{A}^{-1}\)) puis à la multiplier par le vecteur \(\mathbf{b}\) pour obtenir le vecteur solution \(\mathbf{x}\). Mathématiquement, cela peut être exprimé comme suit :

Cette approche nécessite que la matrice \(\mathbf{A}\) soit carrée (c'est-à-dire qu'elle ait le même nombre de lignes et de colonnes) et inversible (c'est-à-dire que son déterminant soit non nul). Le processus implique deux étapes principales :

1. Calculer l'inverse de \(\mathbf{A}\) : \(\mathbf{A}^{-1}\).
2. Multiplier \(\mathbf{A}^{-1}\) par le vecteur \(\mathbf{b}\).

Voici l'équation matricielle et sa solution en détail :

Il est important de noter que cette méthode n'est pas toujours pratique pour les grands systèmes ou les systèmes où \(\mathbf{A}\) n'est pas carrée. En fait, si la taille de ce système est grande, calculer son inverse est coûteux en calcul. Dans de tels cas, d'autres méthodes numériques sont préférées.

Pour calculer l'inverse d'une matrice avec python, l'exemple suivant pourrait être utile :

import numpy as np
from scipy import linalg

a = np.array([[1., 2.], [3., 4.]])
inverse_of_a = linalg.inv(a)

Matrice triangulaire

Une matrice triangulaire est un type spécial de matrice carrée où tous les entrées au-dessus ou en dessous de la diagonale principale sont nulles.

Types de matrices triangulaires :

1. Matrice triangulaire supérieure : Toutes les entrées en dessous de la diagonale principale sont nulles. \(U = \left[\begin{array}{cccc} u_{11} & u_{12} & \cdots & u_{1n} \\ 0 & u_{22} & \cdots & u_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & u_{nn} \end{array}\right]\)

2. Matrice triangulaire inférieure : Toutes les entrées au-dessus de la diagonale principale sont nulles. \(L = \left[\begin{array}{cccc} l_{11} & 0 & \cdots & 0 \\ l_{21} & l_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & \cdots & l_{nn} \end{array}\right]\)

Résoudre un système d'équations linéaires \(A\mathbf{x} = \mathbf{b}\) est considérablement simplifié lorsque \(A\) est une matrice triangulaire car les équations peuvent être résolues séquentiellement.

Pour une matrice triangulaire supérieure \(U\), le système peut être résolu en utilisant la substitution arrière :

Donné : \(U\mathbf{x} = \mathbf{b}\)

1. Commencez par la dernière équation \(u_{nn}x_n = b_n\) pour résoudre \(x_n\): \(x_n = \frac{b_n}{u_{nn}}\)

2. Substituez \(x_n\) dans l'avant-dernière équation pour résoudre \(x_{n-1}\): \(u_{n-1,n-1}x_{n-1} + u_{n-1,n}x_n = b_{n-1}\) \(x_{n-1} = \frac{b_{n-1} - u_{n-1,n}x_n}{u_{n-1,n-1}}\)

3. Continuez ce processus vers le haut pour résoudre \(x_{n-2}, x_{n-3}, \ldots, x_1\).

Ainsi, lorsque \(A\) est une matrice triangulaire, le système \(A\mathbf{x} = \mathbf{b}\) peut être résolu efficacement par substitution arrière (pour les matrices triangulaires supérieures) ou substitution avant (pour les matrices triangulaires inférieures).

Si votre matrice est triangulaire, vous pouvez calculer la solution du système linéaire en utilisant solve_triangular :

import numpy as np
from scipy.linalg import solve_triangular

a = np.array([[3, 0, 0, 0], [2, 1, 0, 0], [1, 0, 1, 0], [1, 1, 1, 1]])
b = np.array([4, 2, 4, 2])

x = solve_triangular(a, b, lower=True)

Élimination de Gauss

L'élimination de Gauss est une méthode utilisée pour résoudre les systèmes d'équations linéaires. Elle transforme un système en une forme de matrice triangulaire supérieure, qui peut ensuite être facilement résolue par substitution arrière. Voici une description étape par étape du fonctionnement de l'algorithme d'élimination de Gauss :

1. Représenter le système comme une matrice augmentée :

   Considérez un système d'équations linéaires :

   \[\begin{cases} a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n = b_1 \\ a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n = b_2 \\ \vdots \\ a_{m1}x_1 + a_{m2}x_2 + \cdots + a_{mn}x_n = b_m \\ \end{cases}\]

   Ceci peut être écrit comme une matrice augmentée :

   \[\left[\begin{array}{cccc|c} a_{11} & a_{12} & \cdots & a_{1n} & b_1 \\ a_{21} & a_{22} & \cdots & a_{2n} & b_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} & b_m \end{array}\right]\]

2. Élimination directe

   L'objectif de l'élimination directe est de transformer la matrice augmentée en une forme triangulaire supérieure.

   Étape 2.1 : Pivotement partiel (optionnel mais recommandé) Pour améliorer la stabilité numérique, échangez les lignes de manière à ce que l'élément le plus grand (en valeur absolue) dans la colonne soit placé dans la position de pivot.

   Étape 2.2 : Éliminer les entrées sous le pivot :

   Pour chaque ligne \(i\) (à partir de la première ligne) :

   • Identifiez l'élément pivot, \(a_{ii}\).
   • Pour chaque ligne subséquente \(j > i\) :
     • Calculez le facteur \(f = \frac{a_{ji}}{a_{ii}}\).
     • Mettez à jour la ligne \(j\) : \(\text{Ligne}_j \leftarrow \text{Ligne}_j - f \times \text{Ligne}_i\).

   Ce processus annule les éléments sous le pivot.

3. Substitution arrière.

   Une fois que la matrice est sous forme triangulaire supérieure, la substitution arrière est utilisée pour trouver les valeurs des variables.

Exemple

Considérez le système d'équations :

1. Matrice augmentée:

   \[\left[\begin{array}{cc|c} 2 & 3 & 8 \\ 4 & 5 & 14 \end{array}\right]\]
2. Élimination directe:
   • Élément pivot : 2 (première ligne, première colonne).
   • Éliminez l'élément dans la deuxième ligne, première colonne : \(f = \frac{4}{2} = 2\). Mettez à jour la deuxième ligne : \(\text{Ligne}_2 = \text{Ligne}_2 - 2 \times \text{Ligne}_1 = [4, 5, 14] - 2 \times [2, 3, 8] = [0, -1, -2]\)

   La matrice devient : \(\left[\begin{array}{cc|c} 2 & 3 & 8 \\ 0 & -1 & -2 \end{array}\right]\)

3. Substitution arrière:
   • À partir de la deuxième ligne : \(-y = -2 \Rightarrow y = 2\).
   • Substituez \(y = 2\) dans la première ligne : \(2x + 3(2) = 8 \Rightarrow 2x + 6 = 8 \Rightarrow 2x = 2 \Rightarrow x = 1\)

La solution est \(x = 1\), \(y = 2\).

L'élimination de Gauss a une complexité temporelle de \(O(n^3)\), où \(n\) est le nombre d'équations. Pour les très grands systèmes, cela peut être coûteux en calcul.

Décomposition LU

La décomposition LU est une méthode de décomposition d'une matrice \(A\) en produit d'une matrice triangulaire inférieure \(L\) et d'une matrice triangulaire supérieure \(U\).

Étant donné une matrice \(A\), décomposez-la en \(A = LU\), où :

• \(L\) est une matrice triangulaire inférieure avec des uns sur la diagonale.
• \(U\) est une matrice triangulaire supérieure.

Une fois décomposée, la solution du système linéaire peut être calculée comme suit :

Pour \(Ax = b\):

1. Décomposez \(A\) en \(LU\).
2. Résolvez \(Ly = b\) pour \(y\) en utilisant la substitution avant.
3. Résolvez \(Ux = y\) pour \(x\) en utilisant la substitution arrière.

Décomposition de Cholesky

La décomposition de Cholesky d'une matrice hermitienne positive définie A, est une décomposition de la forme

où L est une matrice triangulaire inférieure avec des entrées diagonales réelles et positives, et \(L^*\) désigne la transposée conjuguée de L. Chaque matrice hermitienne positive définie (et donc aussi chaque matrice symétrique positive définie réelle) a une décomposition de Cholesky unique.

Une variante étroitement liée à la décomposition classique de Cholesky est la décomposition LDL,

où L est une matrice triangulaire inférieure unitaire (unitriangulaire), et D est une matrice diagonale. C'est-à-dire, les éléments diagonaux de L sont requis pour être 1 au prix de l'introduction d'une matrice diagonale supplémentaire D dans la décomposition.

Si A est symétrique et positive définie, alors \(A x = b\) peut être résolu en calculant d'abord la décomposition de Cholesky \(A = L \times L^*\), puis en résolvant \(L y = b\) pour y par substitution avant, et enfin en résolvant \(L x = y\) pour x par substitution arrière.

Pour calculer la décomposition de Cholesky d'une hermitienne positive définie, utilisez la fonction cholesky de scipy.linalg :

import numpy as np
from scipy.linalg import cholesky

a = np.array([[1,-2j],[2j,5]])
L = cholesky(a, lower=True)

Pour les systèmes linéaires qui peuvent être mis sous forme symétrique, la décomposition de Cholesky (ou sa variante LDL) est la méthode de choix, pour une efficacité supérieure et une stabilité numérique. Comparée à la décomposition LU, elle est environ deux fois plus efficace.

Décomposition QR

Si A n'est pas carrée, nous utilisons la décomposition QR. C’est une décomposition générale qui peut être utilisée pour les matrices carrées et non carrées. La décomposition QR, également connue sous le nom de factorisation QR ou factorisation QU, est une décomposition d'une matrice A en produit \(A = QR\) d'une matrice orthonormale Q et d'une matrice triangulaire supérieure R.

Il existe plusieurs méthodes pour calculer réellement la décomposition QR, telles que le processus de Gram-Schmidt ou les transformations de Householder.

Étapes pour résoudre \(Ax = b\) :

1. Factorisation QR de \(A^T\) :

   Tout d'abord, calculez la factorisation QR de la transposition de \(A\) : \(A^T = QR\), où \(Q\) est une matrice orthogonale \(n \times m\) (donc \(Q^T Q = I\)) et \(R\) est une matrice triangulaire supérieure \(m \times m\).

2. Divisez \(R\) :

   \(R\) prend la forme \(R = \begin{bmatrix} R_1 \\ 0 \end{bmatrix}\), où \(R_1\) est une matrice triangulaire supérieure \(m \times m\), et \(0\) est une matrice zéro \((n-m) \times m\) zéro.

3. Résoudre en utilisant la factorisation QR :

   Pour trouver \(x\), utilisez la formule dérivée de la factorisation QR : \(x = Q \begin{bmatrix} (R_1^T)^{-1} b \\ 0 \end{bmatrix}\) Ici, \((R_1^T)^{-1}\) désigne l'inverse de la transposition de \(R_1\).

Nous pouvons soit trouver \(R_1\) par élimination de Gauss, soit calculer \((R_1^T)^{-1} b\) directement par substitution avant.

Méthodes itératives

Les algorithmes itératifs pour résoudre \(Ax = b\) sont employés lorsque des méthodes comme l'élimination de Gauss prennent trop de temps ou nécessitent trop de mémoire. Contrairement aux méthodes directes, les méthodes itératives ne fournissent généralement pas une solution exacte après un nombre fini d'étapes; au lieu de cela, elles réduisent l'erreur de manière incrémentielle à chaque étape.

Le taux de convergence des méthodes itératives est influencé par les propriétés spectrales de la matrice des coefficients. Par conséquent, transformer le système linéaire en un équivalent avec des propriétés spectrales plus favorables peut être bénéfique. Un préconditionneur est une matrice utilisée pour effectuer une telle transformation.

Par exemple, si une matrice approxime la matrice des coefficients d'une certaine manière, le système transformé

a la même solution que le système original \(Ax = b\), mais les propriétés spectrales de la matrice des coefficients \(M^{-1} A\) peuvent être plus favorables. .

Préconditionneur de Jacobi

Le préconditionneur de Jacobi \(M\) est généralement une matrice diagonale construite à partir des éléments diagonaux de \(A\). Plus précisément, \(M=diag(A)\), où \(diag(A)\) désigne les éléments diagonaux de \(A\).

Méthode de Gauss-Seidel

La méthode de Gauss-Seidel pour résoudre \(n\) équations linéaires commence par une approximation initiale \(X_0\) à la solution, qui n'a pas besoin d'être très précise. La méthode affine cette approximation de manière itérative jusqu'à convergence, où la convergence est déterminée par \(\| x^{(m+1)} - x^{(m)} \|_\infty\) étant inférieure à \(\epsilon\).

Algorithme :

1. Définissez \(x = X_0\)
2. Répéter
   • Définissez \(maxdiff = 0\)
   • Pour \(i = 1\) à \(n\) faire
     • Calculez \(y = (b_i - \sum_{\substack{j=1 \\ j \neq i}}^{n} a_{ij} x_j ) / a_{ii}\)
     • Si \(\|y - x_i\| > maxdiff\) alors
       • Définissez \(maxdiff = \|y - x_i\|\)
     • Mettez à jour \(x_i = y\)
   • Fin pour
3. Jusqu'à \(maxdiff < \epsilon\)

\(x\) est la solution affinée.