None .. An html document created by ipypublish
outline: ipypublish.templates.outline_schemas/rst_outline.rst.j2 with segments: - nbsphinx-ipypublish-content: ipypublish sphinx content
2. Systèmes Linéaires¶
Basile Marchand — Mines ParisTech
2.1. Résolution de systèmes linéaires : Quel intérêt ?¶
2.1.1. Les systèmes linéaires : incontournables¶
En mécanique numérique (non-)linéaire la résolution de systèmes linéaires est impérative
Pour un comportement linéaire :
où \(\mathbf{K}\in\mathbb{R}^{n_{dof}\times n_{dof}}\)
Pour un comportement non-linéaire :
En introduisant un algorithme de Newton ce problème se ramène à :
avec
2.1.2. Notation pour la suite¶
Dans toute la suite du cours on cherchera à résoudre le système linéaire :
Avec : * La matrice \(\mathbf{A}\in \mathbb{R}^{n \times n}\) * Le second membre \(\mathbf{b} \in\mathbb{R}^{n}\) * Le vecteur d’inconnues \(\mathbf{x} \in \mathbb{R}^{n}\)
2.1.3. Les deux classes de méthodes¶
Pour résoudre un système linéaire
Il existe deux grandes classes de méthodes : 1. Les méthodes dites directes 2. Les méthodes dites itératives
Quelle différence ?
Avec les méthodes directs on calcul une solution exacte tandis qu’avec des méthodes itératives on va construire une solution approchée
2.2. Quelques méthodes directes¶
2.2.1. Inverse et factorisation¶
On note souvent :
Mais dans la pratique on ne calcul jamais la matrice inverse \(\mathbf{A}^{-1}\).
Alors comment on fait ?
On calcul une décomposition de \(A\) permettant une résolution rapide
Pour cela on se ramène à des matrices triangulaires
2.2.2. Les matrices triangulaires¶
[1]:
def resolution_triangulaire_superieure(matrice, vecteur):
sol = np.zeros_like(vecteur)
n = vecteur.shape[0]
for i in range(n-1,-1,-1):
sol[i] = vecteur[i]
for j in range(i+1,n):
sol[i] -= matrice[i,j]*sol[j]
sol[i] /= matrice[i,i]
return sol
2.2.3. Élimination de Gauss¶
Règles :
On peut permuter deux lignes de \(\mathbf{A}\) sans changer la solution \(\mathbf{x}\), en effectuant la même permutation de ligne sur \(\mathbf{B}\).
On peut modifier une ligne de \(\mathbf{A}\) avec une combinaison linéaire d’autres lignes de \(\mathbf{A}\) sans changer la solution \(\mathbf{x}\), en effectuant la même com binaison de lignes sur \(\mathbf{B}\)
On peut permuter deux colonnes de \(\mathbf{A}\) si l’on permute les lignes correspondantes de \(\mathbf{x}\). Attention, il faut reconstruire la solution avec la numérotati on d’origine.
On va donc transformer le système
en un autre système
où \(\mathbf{\tilde{A}}\in\mathbb{R}^N\) est une matrice triangulaire supérieure
[2]:
def pivot_gauss(lhs, rhs):
n = rhs.shape[0]
for i in range(0,n):
for j in range(i+1,n):
coef = lhs[j,i]/lhs[i,i]
lhs[j,:] -= coef * lhs[i,:]
rhs[j,0] -= coef * rhs[i,0]
Défauts :
Approche coûteuse
Ne peux servir que pour un seul second membre, i.e. n’est pas réutilisable
2.2.4. La décomposition LU (Lower Upper)¶
Quel intérêt ?
On introduit alors \(\mathbf{y} = \mathbf{U}\cdot \mathbf{x}\)
On calcul alors \(\mathbf{y}\) par une descente :
Et on en déduit \(\mathbf{x}\) par une remontée :
La décomposition LU n’est pas unique, il faut fixer une des deux diagonales à l’identité.
On procède par identification du produit \(\mathbf{L}\cdot \mathbf{U} = \mathbf{A}\)
Cela se traduit par :
Concrètement, pour une matrice \(\mathbf{A}\in\mathbb{R}^{n\times n}\) la factorisation LU se détermine de la manière suivante :
[3]:
def decomposition_lu( lhs ):
L = np.zeros_like( lhs )
U = np.zeros_like( lhs )
n = lhs.shape[0]
for j in range(n,-1,-1):
U[j,j] = 1.
for i in range(j,n):
L[i,j] = A[i,j]
for k in range(j+1):
L[i,j] -= L[i,k]*U[k,j]
for i in range(j+1,n):
tmp = A[j,i]
for k in range(j+1):
tmp -= L[j,k]*U[k,i]
if L[j,j] < 1.e-10:
## error
return
U[j,i] = tmp / L[j,j]
return L, U
Avantages :
Moins coûteux que l’approche précédente
La factorisée est réutilisable
2.2.5. La factorisation de Cholesky¶
On calcul alors la solution de \(\mathbf{A}\cdot \mathbf{x} = \mathbf{b}\) de la manière suivante :
\(\mathbf{L} \cdot \mathbf{y} = \mathbf{b} \; \Rightarrow \; \mathbf{y}\)
\(\mathbf{L}^T \cdot \mathbf{x} = \mathbf{y} \; \Rightarrow \; \mathbf{x}\)
Quel intérêt par rapport à un LU ?
Calcul de la décomposition moins coûteux
On ne stocke que \(\mathbf{L}\) en mémoire
[4]:
def decomposition_cholesky( lhs ):
L = np.zeros_like( lhs )
n = lhs.shape[0]
for i in range(n):
for j in range(i+1):
tmp =0
for k in range(i+1):
tmp += L[i,k] * L[j,k]
if j==i:
L[i,j] = (lhs[i,i] - tmp)**(1./2.)
else:
L[i,j] = 1./L[j,j] * (rhs[i,j] - tmp)
return L
2.2.6. Dans le cadre Python¶
import numpy as np
import numpy.linalg as npl
Deux fonctions utiles :
npl.solve(lhs, rhs)npl.cholesky( lhs )
Dans scipy.linalg un peu plus de choses disponibles :
import scipy.linalg as scl
scl.lu(lhs, rhs): résolution par LUscl.lu_factor( lhs ): décomposition LU
2.3. Les méthodes itératives¶
2.3.1. Alternative aux méthodes directes¶
Comme dit en introduction il y a un second type de méthode pour résoudre le système
Il s’agit des méthodes itératives
Avantage
Moins coûteuses en terme d’occupation mémoire
Inconvénient
Convergence plus ou moins longue dépendant de beaucoup de paramètres méthode employée, conditionnement de la matrice, initialisation
2.3.2. Méthode de Gauss-Seidel¶
L’idée est d’écrire la matrice \(\mathbf{A}\) comme la somme de trois matrices :
Avec \(\mathbf{L}\) triangulaire inférieure, \(\mathbf{U}\) triangulaire supérieure et \(\mathbf{L}\) diagonale. On peut alors construire la suite :
Cet algorithme est associé au problème :
C’est un problème de point fixe. Il y a donc des conditions nécessaires de convergence !
2.3.3. Méthodes de Krylov¶
Deux grandes méthodes majoritairement employées :
Le gradient conjugué
Uniquement pour les matrices symétriques
Le GMRes
Fonctionne pour le non-symétrique
Le principe construire des directions de descentes successive qui soient \(\mathbf{A}\) orthogonale, permet d’améliorer la convergence.
Dans la pratique les solveurs itératifs offres beaucoup d’avantages :
Coût de calcul réduit
Occupation mémoire réduite
Facielement parallélisable
Ils sont en revanche un défaut majeur c’est que la convergence du solveur est directement liée au conditionnement du système linéaire. Or dans la cadre de calcul éléments finis non-linéaires le système linéarisé est très généralement mal conditionné.
Il existe néanmoins des méthodes permettant d’améliorer la convergence des solveurs itératifs, cela passe par la mise en place d’un pré-conditionneur.
Les préconditionneurs classiques sont :
\(\text{diag}( \mathbf{A} ) ^{-1}\)
Factorisation par bloques
Factorisation incomplète
2.3.4. Dans Python¶
Dans les modules de matrices dense aucun solveurs itératif n’est implémenté. En revanche un grand nombre sont disponibles dans `scipy.sparse.linalg’, nous reviendrons sur ce module un peu plus tard.
2.4. Les arrondis numériques¶
2.4.1. Le conditionnement d’une matrice¶
Par définition :
Considérons l’erreur \(\delta \mathbf{A}\) sur la matrice:
avec \(\delta \mathbf{A} \cdot \delta \mathbf{x}\) négligeable devant \(\mathbf{A}\cdot \delta \mathbf{x}\). Ainsi:
Donc :
En conclusion:
Un mauvais conditionnement amplifie les erreurs d’arrondi générées lors de la décomposition de \(\mathbf{A}\).
Exemple d’une matrice mal conditionnée
Observons l’impact d’une perturbation de \(\mathbf{b}\)
[5]:
import numpy as np
A = np.array([[ 7, 1, 11, 10.],[2.,6,5,2.],[8,11,3,8.],[6,9,3,6]])
b = np.array([29.,15,30,24]).reshape((-1,1))
x = np.linalg.solve( A, b)
print(" cond( A ) = {}".format( np.linalg.cond(A) ) )
cond( A ) = 1424.9502711293317
[6]:
delta = [0.01, 0.03, 0.05, 0.07, 0.09, 0.1, 0.3, 0.5, 0.7, 0.9, 1.]
perturb = np.array([1.,-1.,1.,-1.]).reshape((-1,1))
error_b = [0.]*len(delta)
error_x = [0.]*len(delta)
norm_b = np.linalg.norm(b)
norm_x = np.linalg.norm(x)
for i, d in enumerate(delta):
error_b[i] = np.linalg.norm( d * perturb ) / norm_b
xi = np.linalg.solve(A, b+d*perturb )
error_x[i] = np.linalg.norm( x - xi ) / norm_x
import matplotlib.pyplot as plt
%matplotlib inline
plt.plot( error_b, error_x )
[6]:
[<matplotlib.lines.Line2D at 0x7f9617b80350>]
2.4.2. Améliorer le conditionnement¶
On peut par exemple définir la matrice de scaling comme :
2.5. Quelques mots sur la compléxité algorithmique¶
2.5.1. La compléxité d’un algorithme¶
Problématique : Comment évaluer les performances d’un algorithme ? Solution 1 : Je l’implémente et je mesure le temps Problème : Les performances dépendent alors de la configuration de l’ordinateur
Si je donne à mon algorithme une entrée de taille N, quel est l’ordre de grandeur en fonction de N du nombre d’opération que je vais effectuer ?”
On note la compléxité \(T(N)\)
2.5.2. Exemple de complexité¶
Faire la somme de tous les termes d’une liste
Faire une moyenne pondérée
Trier une liste :
Tri à bulle : \(T(N) = O(N^2)\)
Tri rapide : \(T(N) = O(N \log N)\)
Produit Matrice-Vecteur
2.5.3. Pour la résolution des systèmes linéaires¶
Résolution d’un système triangulaire :
Pivot de Gauss :
Décomposition LU :
Décomposition de Cholesky :
\(n^3/6\) additions+multiplications, \(n(n-1)/2\) divisions, \(n\) racines carrés
2.6. Le cas particulier de nos matrices¶
Il ne faut pas oublier que les systèmes que l’on cherche à résoudre ne sont pas quelconques ils proviennent d’une formulation éléments finis. Qu’est ce que cela implique me répondrez vous. Beaucoup de choses… Tout d’abord le fait que l’on utilise une discrétisation éléments finis implique que nos matrices sont nécéssairement creuses, i.e. elles comprennent un grand nombre de valeurs nulles. Cela vient du fait que les fonctions de forme utilisées dans la méthode des éléments finis sont à support compact.
2.6.1. Les matrices sparses¶
Ne stocker que les valeurs non nulles de la matrices
Différents formats : COO, CSR, CSC
On peut alors tirer parti de ces formats de stockage pour réduire la complexité des algorithmes d’algèbre linéaire
Par exemple un produit matrice-vecteur en sparse est de complexité :
Avec \(\omega\) la largeur de bande de la matrice
2.7. Exercice — Système masses-ressorts¶
[7]:
from IPython.display import IFrame
IFrame("media/spring/spring.pdf", width=600, height=300)
[7]:
<IPython.lib.display.IFrame at 0x7f9617e4c590>
L’objectif est de déterminer par une approche matricielle la résponse du système à un effort \(F\).
La formulation du problème doit donc se ramener à la résolution d’un problème de la forme
Avec \(\mathbf{K} \in \mathbb{R}^{M\times M}\), \(\mathbf{u} \in \mathbb{R}^{M}\), \(\mathbf{F} \in \mathbb{R}^{M}\) et \(M\) le nombre de points du système.
Pour cela on rappel que l’énergie potentielle du système peut s’exprimer de la forme suivante :
Avec \(u_{i,0}\) le déplacement du premier noeuds d’attache du \(i\)-ème ressort et \(u_{i,1}\) le déplacement du second noeud d’attache du \(i\)-ème ressort.
Cette énergie potentielle peut s’écrire sous la forme matricielle suivante :
En utilisant alors le théorème de l’énergie potentielle nous pouvons écrire que
Question 1 :
A partir de l’expression de l’énergie potentielle d’un ressort définir la matrice de rigidité élémentaire.
Question 2 :
Utiliser la matrice de rigidité d’un ressort pour construire la matrice \(\mathbf{K}\) globale. Pour cela utiliser la notion de table de connectivité. Pour rappel la table de connectivité est une liste \(L\) telle que le \(i\)-ème élément de \(L\) est le doublet des indices des points d’attache du ressort.
Question 3 :
Construire le second membre \(F\).
Question 4 :
Résoudre le système linéaire
Bonus :
En utilisant matplotlib, visualiser le profil de la matrice \(\mathbf{K}\):
import matplotlib.pyplot as plt
plt.imshow( K )
plt.colorbar()
plt.show()
Que peut-on en conclure ? Quelle piste d’amélioration est envisageable pour accélérer la résolution du problème ?
[ ]: