On montre, pour commencer, comment faire pour résoudre numériquement un système linéaire, en utilisant numpy.
On commence par définir la matrice $A$ et le vecteur $b$ du système $AX=b$ qui nous intéresse :
>>> from numpy import * >>> A = matrix([[1,2,3],[4,5,6],[7,8,9]]) >>> A matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) >>> b=matrix((1,2,3)).reshape((3,1)) >>> b matrix([[1], [2], [3]])
On importe ensuite le paquetage numpy.linalg, qui contient la fonction solve permettant cette résolution, et on calcule $X$ :
>>> from numpy.linalg import * >>> X=solve(A,b) >>> X matrix([[-0.26193963], [ 0.52387926], [ 0.0713937 ]])
Pour finir, on vérifie que le produit $AX$ redonne bien $b$...
>>> A*X matrix([[ 1.], [ 2.], [ 3.]])
Pour résoudre le système $Ax = b$, où $A$ est une matrice quelconque, on cherche à se ramener à un système plus simple que l'on sait résoudre : un système triangulaire (méthode du pivot de Gauss) ou un système diagonal (méthode de Gauss-Jordan) Pour cela on transforme un système en un système équivalent obtenu par :
Soit à résoudre le système suivant de trois équations à trois inconnues
$ \left\{ \begin{array}{rrrrrrrr} x & +& 2y& + & 2z& = & 2 &L_1\\ x & +& 3y& - & 2 z & = & -1&L_2\\ 3x & +& 5y & + & 8z & = & 8 &L_3\end{array} \right.$
Le premier pivot est $a_{11} = 1 $. On supprime l'inconnue $x$ de $L_2$ et $L_3$ en utilisant la transformation :
On a :
$ \left\{ \begin{array}{rrrrrrrrr} x & +& 2y& + & 2z &= & 2 & & &L'_1\\ & & y& - & 4z &= &-3 & & & L'_2 \leftarrow L_2 - \frac{a_{21}}{a_{11}} \cdot L_1\\ & & -y& + & 2z &= &2 & & & L'_3 \leftarrow L_3 - \frac{a_{31}}{a_{11}} \cdot L_1\end{array} \right.$
Le second pivot est $a'_{22} = 1 $.
On supprime l'inconnue $y$ de $L_3$ en utilisant la transformation :
Ce qui donne, cette fois-ci...
$ \left\{ \begin{array}{rrrrrrrrrr} x & +& 2y& + & 2z &= & 2 & & & L"_1\\ & & y& - & 4z &= &-3 & & & L"_2 \\ & & & & -2z &= &-1 & & & L"_3 = L'_3 - \frac{a'_{32}}{a'_{22}} \cdot L'_2 \end{array} \right.$
Une fois la valeur de $z$ obtenue, il reste à remonter l'information pour obtenir les valeurs des variables $y$ puis $x$.
La méthode de Gauss est itérative : on considère les suites $(A^n)$ et $(b^n)$ de matrices telle que $A^0 = A$ et $b^0 = b $ (attention à la notation utilisée : $X^n$ n'est pas une nième puissance, mais le nième terme d'une suite).
On construit successivement la séquence de matrices à $n$ ligne et $n+1$ colonnes $(A^0, b^0)$, $(A^1, b^1)$, ..., $(A^n, b^n)$ telle qu'à l'itération $i$, la matrice $(A^i,b^i)$ a la forme
$\left(\begin{array}{rrrr|rrrr|r} a_{1,1}^i & a_{1,2}^i & \ldots & a_{1,i}^i & a_{1,i+1}^i & \ldots & a_{1,n-1}^i & a_{1,n}^i & b_1^i\\ 0 & a_{2,2}^i & \ldots & a_{2,i}^i & a_{2,i+1}^i & \ldots & a_{2,n-1}^i & a_{2,n}^i &b_2^i \\ \vdots & \ddots & \ddots & \vdots & & \vdots & \vdots & \vdots & \vdots \\ 0 & \ldots & 0 & a_{i,i}^i & a_{i,i+1}^i & \ldots & a_{i,n-1}^i & a_{i,n}^i & b_n^i \\0 & 0 & \ldots & O & a_{i+1,i+1}^i& \ldots & a_{i+1,n-1}^i & a_{i+1,n}^i & b_{n+1}^i\\ \vdots & \vdots & & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \ldots & O & a_{n-1,i+1}^i& \ldots & a_{n-1,n-1}^i & a_{n-1,n}^i & b_{n-1}^i\\ 0 & 0 & \ldots & O & a_{n,i+1}^i & \ldots & a_{n,n-1}^i & a_{n,n}^i & b_{n}^i\\ \end{array}\right)$
On obtient la matrice $A^i$ en faisant apparaître des 0 à la colonne $i+1$. Pour cela chaque ligne $L_j$, $j > i+1$, est remplacée par une combinaison linéaire d'elle même avec la ligne $L_{i+1}$ comme suit :
$L_j = L_j - \frac{a_{i+1,j}^i}{a_{i+1,i+1}^i}\cdot L_{i+1}$ si $a_{i+1,i+1}^i$ est non nul.
C'est une variante de la méthode du pivot de Gauss : à la l'étape $i$, on combine toutes les lignes (sauf la ligne i) avec la ligne $i$ (au lieu de ne le faire que pour les lignes d'indice supérieurs à $i$). L'intéret est que cela fait apparaître des 0 sur toute la colonne sauf au niveau du pivot $a_{ii}^i$
Soit à résoudre le système de trois équations à trois inconnues suivant :
$ \left\{ \begin{array}{rrrrrrrl} 2x & +& y& - & 4z& = & 8 &L_1\\ 3x & +& 3y& - & 5z & = & 14&L_2\\ 4x & +& 5y & - & 2z & = & 16 &L_3\end{array} \right.$
On trouve alors...
$ \left\{ \begin{array}{rrrrrrrl} x & +& 1/2y& - & 2z &= & 4 & L'_1 \leftarrow L_1/2\\ & & 3/2 \cdot y& + & z &= &2 & L'_2 \leftarrow L_2 -3 \cdot L'_1\\ & & 3y& + & 6z &= &0 & L'_3 \leftarrow L_3 -4 \cdot L'_1 \end{array}\right.$
Cela donne :
$ \left\{ \begin{array}{rrrrrrrl} x & & & & -7/3 \cdot z &= & 10/3 & L"_1 \leftarrow L'_1 -1/2 \cdot L"_2 \\ & & y& + & 4/3 \cdot z &= &4/3 & L"_2 \leftarrow 2/3 \cdot L'_2 \\ & & & & 4z &= &-4 & L"_3 \leftarrow L'_3 -3 \cdot L"_2 \end{array}\right.$
D'où, finalement,
$ \left\{ \begin{array}{rrrrrrrl} x & & & & & = & 1 & L_1^{3} \leftarrow L"_1 + 7/3 \cdot L_2^{3}\\ & & y& & &= &2 & L_2^{3} \leftarrow L"_2 - 2/3 \cdot L_2^{3} \\ & & & & z &= &-1 & L_3^{3} \leftarrow L"_3/4 \cdot L"_2 \end{array}\right.$
et on trouve directement les racines sur la 4ème colonne.
On a le pseudo algorithme suivant :
pour i = 1 à n recherche du pivot (non nul) sur la ligne p, p >= i échange éventuel des lignes i et p -> le pivot est Aii division de la ligne Li par Aii pour obtenir L'i pour k = 1 à n sauf i L'k = L_k - Aki .L'i
Les solutions sont dans la colonne $n+1$.
$\left(\begin{array}{cc}10^{-12} & 1 \\ 1 & 1 \end{array}\right) \cdot \left(\begin{array}{c}x \\y \end{array}\right)=\left(\begin{array}{c} 1 \\2 \end{array}\right)$ et $\left(\begin{array}{cc}1 & 1 \\ 10^{-12} & 1 \end{array}\right) \cdot \left(\begin{array}{c}x \\y \end{array}\right)=\left(\begin{array}{c} 1 \\2 \end{array}\right)$.
$ \left(\begin{array}{cccc}10 & 7 & 8 & 7 \\7 & 5 & 6 & 5 \\9 & 6 & 10& 9 \\ 7 & 5 & 9 & 10 \end{array} \right)$
$\left( \begin{array}{c} 32 \\ 23 \\ 33 \\ 31 \end{array} \right)$ $\left( \begin{array}{c} 32,1 \\ 22,9 \\ 33,1 \\ 30,9 \end{array} \right)$.