CAMA : ma10 Système d'équations

Cours du 26 / 04

Systèmes matriciels

Un système de plusieurs équations à autant d'inconnues peut s'écrire comme système matriciel

$A\mathbf{x}=\mathbf{b}$ :

$$\left[\begin{array}{c}{a}_{11}{a}_{12}\dots a_{1n}\\ a_{21}{a}_{22}\dots a_{2n}\\ ⋮\\ a_{n1}{a}_{n2}\dots a_{nn}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right]$$

Exemple

On acheté 3 fois des quantités de fruits dont nous n'avons pas le prix.

# A est la quantité de chaque fruit achetée
# x est le prix de chaque fruit
# b est la somme qu'on a payé pour chaque course

A = np.array([[6,5,4], [5,3,2], [7,3,2]])
b = np.array([11.7, 7.9, 9.5])

x = lin.inv(A) @ b

array([0.8, 0.9, 0.6])

Résolution d'un système matriciel

Méthode du pivot de Gauss

On transforme

$A$ en une matrice triangulaire pour résoudre le système de

$O(n^2)$ opérations.

On met des

$0$ sur la première colonne en dessous de la diagonale en multipliant

$A$ par la matrice

$E_1$ suivante :

$$E_1=\left[\begin{array}{c}100\dots 0\\ \frac{-a_{21}}{a_{11}}10\dots 0\\ \frac{-a_{31}}{a_{11}}01\dots 0\\ ⋮\\ \frac{-a_{n1}}{a_{11}}00\dots 1\end{array}\right]$$

$E_2$ sera similaire avec des termes

$\frac{-a_{k2}}{a_{22}}$ sous la diagonale, de même pour les matrices

$E_n$ suivantes.

Si on multiplie

$A$ par

$E_1$ il faut également multiplier

$b$ par

$E_1$.

Système matriciel avec matrice triangulaire

Il faut résoudre

$U\mathbf{x}=c$ avec

$U$ une matrice triangulaire supérieure.

On part de la dernière ligne et on obtient

$$\mathbf{x}[-1]=\frac{c[-1]}{U[-1,-1]}$$
Une fois

$\mathbf{x}[-1]$ connu, on en déduit la valeur de

$\mathbf{x}[-2]$, puis celle de

$\mathbf{x}[-3]$, etc.

def solve_gauss(A, b):
    for i in range(len(A)-1):
        E = np.diag(np.array([1.,] * len(A), dtype=A.dtype))
        coefs = - A[i+1:,i] / A[i,i]
        E[i+1:,i] = coefs
        A[i:, i:] = E[i:,i:] @ A[i:,i:]
        b[i+1:] += coefs * b[i]   # multiplication terme à terme
    # A est maintenant triangulaire supérieur
    res = np.zeros(len(b), dtype=b.dtype)
    res[-1] = b[-1] / A[-1,-1]
    for i in range(len(A)-1)[::-1]:
        res[i] = (b[i] - A[i,i+1:] @ res[i+1:]) / A[i,i]
    return res

A = 10 * np.random.random(size=(5,5))
b = A.sum(axis=1)

A
 [[5.655 3.042 3.18  9.672 8.761]
 [3.963 9.923 9.868 7.934 6.328]
 [0.697 9.189 7.799 2.046 2.184]
 [6.314 8.533 4.879 8.112 4.583]
 [3.803 6.89  2.266 6.087 0.361]] 
b
 [30.31  38.015 21.914 32.42  19.407]

solve_gauss(A, b)

Décomposition Lower Upper

Si on a besoin de résoudre plusieurs systèmes matriciels avec

$A$, on décompose

$A$ en un produit d'une matrice triangulaire inférieure et d'une matrice triangulaire supérieure.

$$A=LU$$

On utilise le pivot de Gauss mais au lieu de modifier

$b$, on calcule l'inverse de la matrice

$E_n\dots E_2{E}_1$ :

$$\begin{array}{rl}{E}_n\dots E_2{E}_{1}Ax& =E_n\dots E_2{E}_{1}b\text{donc}\\ (E_n\dots E_2{E}_1{)}^{-1}{E}_n\dots E_2{E}_{1}Ax& =b\\ E_1^{-1}{E}_2^{-1}\dots E_n^{-1}{E}_n\dots E_2{E}_{1}Ax& =b\end{array}$$

Ce calcul est simple, la matrice inverse a les valeurs opposées :

$$E_1^{-1}=\left[\begin{array}{c}100\dots 0\\ \frac{a_{21}}{a_{11}}10\dots 0\\ \frac{a_{31}}{a_{11}}01\dots 0\\ ⋮\\ \frac{a_{n1}}{a_{11}}00\dots 1\end{array}\right]$$

Le produit

$E^{-1}=E_1^{-1}{E}_1^{-1}{E}_2^{-1}{E}_3^{-1}\dots E_n^{-1}$ est la concaténation des colonnes :

$$E^{-1}=\left[\begin{array}{c}100\dots 0\\ \frac{a_{21}}{a_{11}}10\dots 0\\ \frac{a_{31}}{a_{11}}\frac{a_{32}}{a_{22}}1\dots 0\\ ⋮\\ \frac{a_{n1}}{a_{11}}\frac{a_{n2}}{a_{22}}\frac{a_{n3}}{a_{33}}\dots 1\end{array}\right]$$

Pour résoudre

$LU\mathbf{x}=\mathbf{b}$ :

1. on résoud
   $L\mathbf{y}=\mathbf{b}$ obtenant
   $\mathbf{y}$
2. on résoud
   $U\mathbf{x}=\mathbf{y}$ obtenant la solution
   $\mathbf{x}$

def LU(A):
    L = np.diag([1.,] * len(A))
    for i in range(len(A)-1):
        E = np.diag([1.,] * len(A))
        E[i+1:,i] = - A[i+1:,i] / A[i,i]
        L[i+1:,i] = -E[i+1:,i]
        A[i:, i:] = E[i:,i:] @ A[i:,i:]
    return L, A

A = 10 * np.random.random(size=(5,5))
L,U = LU(A.copy())  # Attention, notre fonction modifie A donc si on veut le réutiliser il faut une copie

A
 [[2.697 6.265 5.876 1.927 3.951]
 [2.495 0.021 9.085 0.504 9.23 ]
 [0.788 1.982 5.048 9.656 8.581]
 [3.19  7.474 8.344 0.124 2.577]
 [4.209 6.825 9.223 9.025 4.733]]
L
[[ 1.     0.     0.     0.     0.   ]
 [ 0.925  1.     0.     0.     0.   ]
 [ 0.292 -0.026  1.     0.     0.   ]
 [ 1.183 -0.011  0.418  1.     0.   ]
 [ 1.561  0.511 -0.529 -1.924  1.   ]]
U
[[  2.697   6.265   5.876   1.927   3.951]
 [  0.     -5.776   3.647  -1.279   5.574]
 [  0.      0.      3.427   9.059   7.573]
 [  0.      0.      0.     -5.957  -5.201]
 [  0.      0.      0.      0.    -10.285]]

A - (L @ U)

array([[ 0.,  0.,  0.,  0.,  0.],
       [ 0., -0.,  0.,  0.,  0.],
       [ 0.,  0., -0.,  0.,  0.],
       [ 0.,  0.,  0., -0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.]])

Gauss Jordan

On diagonalise

$A$ par des multiplications matricielles similaire à celles de Gauss qui annulent aussi les termes au dessus de la diagonale :

$$E_3=\left[\begin{array}{c}10\frac{-a_{11}}{a_{33}}0\dots 0\\ 01\frac{-a_{21}}{a_{33}}0\dots 0\\ 0010\dots 0\\ 00\frac{-a_{41}}{a_{33}}1\dots 0\\ ⋮\\ 00\frac{-a_{n1}}{a_{33}}0\dots 1\end{array}\right]$$

def solve_gauss_jordan(A,b):
    for i in range(len(A)):
        d1 = np.diag([1.,] * len(A))
        d1[:,i] = - A[:,i] / A[i,i]
        A = d1 @ A
        b = d1 @ b
    return b / np.diag(A)

A = 10 * np.random.random(size=(5,5))
b = A.sum(axis=1)
solve_gauss_jordan(A, b)

array([1., 1., 1., 1., 1.])