import numpy as np
import numpy.linalg as lin

np.set_printoptions(precision=3, linewidth=150, suppress=True)
np.random.seed(123)

La simulation numérique

Si le prix des bananes est important, le secteur qui a les plus gros besoins en résolution de systèmes matriciels est la simulation numérique. Cela comprend ceci

et tout ce qu'on simule numériquement à savoir tout ce qui touche le transport, l'énergie, la construction, tout ce qu'on fabrique, qui coûte cher et qui doit avoir un comportement physique bien précis. Mais ce n'est pas tout, comprendre notre environnement (météo, réchauffement de la planète, chimie...) revient aussi à résoudre des systèmes matriciels ainsi que d'autres domaines. Cela étant il existe des méthodes de simulation numérique qui ne passent pas par des systèmes matriciels.

Pour faire l'image ci-dessus, on transforme des équations physique comme Navier-Stokes en un système matriciel où les inconnues sont définies en chaque point d'un maillage à définir. Dans notre cas l'inconnue est la pression et le maillage est l'intérieur d'une boite imaginaire qui comprend l'avion et l'air qui circule autour.

Si la boite est un cube avec 1000 points dans chaque direction, on a 1 milliard de points dans la boite (moins ce qui est à l'intérieur de l'avion mais restons sur 1 milliards). Alors la matrice a 1 trillion d'éléments (un milliard au carré).

$$ \begin{bmatrix} a_{11} \; a_{12} \ldots a_{1,10^9} \\ a_{21} \; a_{22} \ldots a_{2,10^9} \\ \vdots \\ a_{10^9,1} a_{n2} \ldots a_{10^9,10^9} \\ \end{bmatrix} \; \begin{bmatrix} p_{1} \\ p_{2} \\ \vdots \\ p_{10^9} \\ \end{bmatrix} = \begin{bmatrix} f_{1} \\ f_{2} \\ \vdots \\ f_{10^9} \\ \end{bmatrix} $$

Inverser une matrice se fait en $O(n³)$ opérations soit $O(10^{27})$ dans notre cas.

Si vous avez l'ordinateur le plus puissant du monde vous pouvez faire 1 Eflops soit $10^{18}$ opérations par seconde en pic (en 2024). Aussi il vous faudra $10^{9}$ secondes soit un peu plus de 30 ans. C'est trop.

Inverser une matrice ou résoudre par une méthode directe n'est pas la bonne solution pour résoudre un grand système matriciel.

Exercice 1

Pour une telle simulation il faut aussi calculer la vitesse de l'air en chacun des 1 milliard de points de notre maillage. Une vitesse c'est 3 variables qu'il faut ajouter à la pression. Combien de temps faut-il pour inverser la matrice qui intègre aussi la vitesse de l'air ?

# ceci est une calculatrice, vous pouvez donc écrire la réponse ici

Estimer le temps d'un gros calcul avant de le lancer est une bonne idée.

D'ailleur avec des temps si longs, il ne faut pas rester aux ordres de grandeurs mais préciser avec la formule exacte. Ainsi avec Choleski c'est $n^3/3$ opérations, donc 3 fois moins, mais bon...

Méthodes itératives

Les méthodes itératives sont des méthodes qui s'approchent pas à pas de la solution recherchée. Elles permettent de trouver une approximation de ${\bf x}$, la solution de $A\, {\bf x} = b$.

En général on arrête le calcul lorsqu'on estime qu'on est à une distance choisie de la solution (qu'on appelle l'erreur) plutôt que d'attendre d'être à la solution exacte. De toute facon la solution exacte sur un ordinateur qui est limité en nombre de chiffres après la virgule peut être inatteignable. Donc on cherchera jamais à avoir une erreur plus petite que notre précision maximale.

L'idée fondamentale de l'algorithme itératif est d'avoir une formule du type $\; {\bf x}^{t+1} = B \, {\bf x}^t + {\bf c}\;$ ou en Python:

x = np.random(size = c.size)
while np.square(x - old_x) > seuil:
    old_x = x
    x = B @ x + c

La magie c'est si x converge. Dans ce cas on a atteint un point fixe c.a.d. que ${\bf x}^{t+1} = {\bf x}^t$ et donc

$${\bf x}^t = B \, {\bf x}^t + {\bf c} \quad \textrm{c.a.d.} \quad (Id -B) \, {\bf x}^t = {\bf c}$$

On retrouve notre $A \; {\bf x} = {\bf b}$ cela étant en pratique on ne découpe pas A en Id et B car ca ne converge pas (sauf cas particuliers). Regardons comment fait la méthode de Jacobi.

Méthode de Jacobi

La méthode de Jacobi découpe la matrice A en M et N avec

• M la matrice diagonale qui comprend les éléments de la diagonale de A
• N = M - A (donc 0 sur la diagonale et l'opposé des éléments de A ailleurs)

Ainsi le système à résoudre est $\; (M - N) \, {\bf x} = {\bf b}$.

La formule itérative est donc pour l'itération k+1 exprimée en fonction de l'itération k :

$$ M \; {\bf x}^{k+1} = N\; {\bf x}^k + {\bf b} $$

et comme M est une matrice diagonale on a :

$$ \begin{array}{l} a_{11} x_1^{k+1} = \qquad -a_{12} \, x_2^k - a_{13} \, x_3^k \ldots - a_{1n} \, x_n^k + b_1 \\ a_{22} x_2^{k+1} = -a_{21} \, x_1^k \qquad - a_{23} \, x_3^k \ldots - a_{2n} \, x_n^k + b_2 \\ a_{33} x_3^{k+1} = -a_{31} \, x_1^k - a_{32} \, x_2^k \qquad \ldots - a_{3n} \, x_n^k + b_3 \\ \vdots \\ a_{nn} x_n^{k+1} = -a_{n1} \, x_1^k - a_{n2} \, x_2^k \ldots - a_{n-1,n-1} \, x_{n-1}^k \qquad + b_n \\ \end{array} $$

On voit en filigramme $A \; {\bf x} = {\bf b}$.

Pour calculer $x_1^{k+1}$ il faut diviser le terme de droite de la première ligne par $a_{11}$ il faut donc que $a_{11} \ne 0$.

Pour calculer les termes suivants $x_i^{k+1}$ il faut aussi que $a_{ii}$ soient non nul donc il faut que A n'ait pas pas de zéro sur sa diagonale.

Cela se retrouve dans l'écriture suivante qui reprend la formule initiale : $ {\bf x}^{k+1} = M^{-1} \;(N\; {\bf x}^k + {\bf b})$

On voit que pour être efficace, il faut que M soit facilement inversible, sinon autant inverser A directement. Ici c'est bien le cas puisque M est diagonale.

Programmons Jacobi

A = np.random.randint(10, size=(4,4))
b = A.sum(axis=1)                     # ainsi la solution est [1,1,1,1]
print('A:\n', A, "\nb:\n", b, "\n")

M = np.diag(A)        # attention, c'est un vecteur
N = np.diag(M) - A    # np.diag d'une matrice donne un vecteur, np.diag d'un vecteur donne une matrice
print(f"M:\n {np.diag(M)}\nN:\n {N}\n")

x = np.random.random(4)
for i in range(20):
    print(f"x_{i} = {x}")
    x = (N @ x + b) / M

A:
 [[2 2 6 1]
 [3 9 6 1]
 [0 1 9 0]
 [0 9 3 4]] 
b:
 [11 19 10 16] 

M:
 [[2 0 0 0]
 [0 9 0 0]
 [0 0 9 0]
 [0 0 0 4]]
N:
 [[ 0 -2 -6 -1]
 [-3  0 -6 -1]
 [ 0 -1  0  0]
 [ 0 -9 -3  0]]

x_0 = [0.398 0.738 0.182 0.175]
x_1 = [4.127 1.837 1.029 2.203]
x_2 = [-0.526 -0.195  0.907 -0.906]
x_3 = [3.427 1.782 1.133 3.759]
x_4 = [-1.56  -0.204  0.913 -0.86 ]
x_5 = [3.395 2.118 1.134 3.775]
x_6 = [-1.907 -0.196  0.876 -1.616]
x_7 = [3.877 2.342 1.133 3.784]
x_8 = [-2.133 -0.357  0.851 -2.12 ]
x_9 = [4.364 2.49  1.151 4.165]
x_10 = [-2.525 -0.574  0.834 -2.467]
x_11 = [4.804 2.671 1.175 4.665]
x_12 = [-3.027 -0.792  0.814 -2.89 ]
x_13 = [5.293 2.898 1.199 5.17 ]
x_14 = [-3.581 -1.027  0.789 -3.421]
x_15 = [5.87  3.159 1.225 5.719]
x_16 = [-4.194 -1.298  0.76  -4.026]
x_17 = [6.531 3.45  1.255 6.35 ]
x_18 = [-4.891 -1.608  0.728 -4.704]
x_19 = [7.277 3.779 1.29  7.073]

Ca ne converge pas vraiment... (si vous relancez et voyez des NaN c'est qu'il y a un zéro sur la diagonale)

2e essai :

A = np.random.randint(10, size=(4,4))
A = A + np.diag(A.sum(axis=0))
b = A.sum(axis=1)                     # ainsi la solution est [1,1,1,1]
print('A:\n', A, "\nb:\n", b, "\n")


M = np.diag(A)        # attention, c'est un vecteur
N = np.diag(M) - A    # np.diag d'une matrice donne un vecteur, np.diag d'un vecteur donne une matrice
print(f"M:\n {np.diag(M)}\nN:\n {N}\n")

x = np.random.random(4)
for i in range(20):
    print(f"x_{i} = {x}")
    x = (N @ x + b) / M

A:
 [[24  2  4  8]
 [ 0 24  9  3]
 [ 4  6 16  5]
 [ 6  2  1 32]] 
b:
 [38 36 31 41] 

M:
 [[24  0  0  0]
 [ 0 24  0  0]
 [ 0  0 16  0]
 [ 0  0  0 32]]
N:
 [[ 0 -2 -4 -8]
 [ 0  0 -9 -3]
 [-4 -6  0 -5]
 [-6 -2 -1  0]]

x_0 = [0.766 0.777 0.028 0.174]
x_1 = [1.456 1.468 1.4   1.088]
x_2 = [0.865 0.839 0.683 0.873]
x_3 = [1.109 1.135 1.134 1.045]
x_4 = [0.951 0.944 0.908 0.967]
x_5 = [1.031 1.039 1.043 1.015]
x_6 = [0.984 0.982 0.973 0.99 ]
x_7 = [1.009 1.011 1.014 1.005]
x_8 = [0.995 0.994 0.992 0.997]
x_9 = [1.003 1.003 1.004 1.002]
x_10 = [0.998 0.998 0.998 0.999]
x_11 = [1.001 1.001 1.001 1.   ]
x_12 = [1.    0.999 0.999 1.   ]
x_13 = [1. 1. 1. 1.]
x_14 = [1. 1. 1. 1.]
x_15 = [1. 1. 1. 1.]
x_16 = [1. 1. 1. 1.]
x_17 = [1. 1. 1. 1.]
x_18 = [1. 1. 1. 1.]
x_19 = [1. 1. 1. 1.]

Magie !

Pourquoi le 2e cas marche ?

Pour qu'une méthode itérative du type ${\bf x} = B\; {\bf x} + {\bf c}$ converge il faut au choix

• $\rho(B) < 1\quad$ avec $\rho$ le rayon spectral (qui est la plus grande valeur propre en valeur absolue)
• $||B|| < 1\quad$ avec une norme matricielle est subordonnée à une norme vectorielle.

Cas de la méthode de Jacobi

Pour la matrice B de Jacobi, on respecte ces conditions si A est à diagonale dominante à savoir que chaque élément de la diagonale est plus grand en module que la sommes tous les autres éléments en module de sa ligne ou colonne ($|a_{i,i}| \ge \sum_{j \ne i} |a_{i,j}|$).

Jacobi converge aussi si A est symétrique, réelle et définie positive (c.a.d. $\forall {\bf x}, \; {\bf x}^T \, A \, {\bf x} > 0$).

Temps calcul

Si on converge en 10 itérations alors on utilise 10 multiplications matricielles, 10 additions vectorielles et 10 divisions vectorielles soit 180 opérations à comparer aux $4^3 / 3 = 22$ opérations d'une méthode directe. Zut !

Quelques raisons qui font qu'on y perd sont

• La matrice A est trop petite, les méthodes itératives sont intéressantes pour les grandes matrice
• La méthode de Jacobi n'est pas la meilleure (mais elle est très facilement parallélisable)
• Le rayon spectral de la matrice est trop grand. Plus le rayon spectral est petit et plus on converge vite.