+++ title = "3. Factorisation LU"
date = 2018-09-09T00:00:00
draft = false # Is this a draft? true/false toc = true # Show table of contents? true/false type = "docs" # Do not modify.
math = true weight = 190 diagram = false #markup = "mmark"
edit_page = {repo_url = "https://github.com/Bertbk/course_4m053", repo_branch = "master", submodule_dir="content/course/4m053/"}
[git] icon = "github" repo = "https://github.com/Bertbk/course_4m053" submodule_dir = "content/course/4m053/"
[menu.4m053] parent = "VI. Solveurs Directs" name = "3. Factorisation LU" weight = 10
+++
- Calculer la factorisation LU d'une matrice
- Résoudre le système linéaire une fois la factorisation effectuée
Cette méthode permet de transformer une matrice carré
Ainsi, comme
- Calculer la factorisation LU
- Résoudre un système linéaire triangulaire inférieur (avec des 1 sur la diagonale)
- Réoudre un système linéaire triangulaire supérieur
C'est parfait, nous avons déjà implémenté [les fonctions de résolution]({{<relref "direct_triang.md">}}), il ne nous manque plus que le calcul de la factorisation LU !
Notons
\begin{equation}
A=\begin{pmatrix}
a_{0,0} & A_{0,1} \\\
A_{1,0} & A_{1,1}
\end{pmatrix} =
\begin{pmatrix}
1 & 0 \\\
L_{1,0} & I
\end{pmatrix}
\begin{pmatrix}
u_{0,0} & U_{0,1} \\\
0 & S_{1,1}
\end{pmatrix}
\label{eq:factorisation_partielle}
\end{equation}
où
{{% alert exercise %}} Vérifiez que :
$L_{1,0} = U_{0,0}^{-1}A_{1,0} = A_{1,0} / u_{0,0}$ $U_{0,1} = A_{0,1}$ -
$S_{1,1}=A_{1,1}-A_{1,0}A_{0,0}^{-1}A_{0,1} = A_{1,1} - L_{1,0}U_{0,1}$ {{% /alert %}}
{{% alert note %}}
Afin d'éviter toute confusion, nous utilisons des lettres et des indices minuscules pour les coefficients (e.g.
{{% alert note %}} La factorisation partielle peut aussi être opérérée par bloc : $$ A=\begin{pmatrix} A_{0,0} & A_{0,1} \\\ A_{1,0} & A_{1,1} \end{pmatrix} = \begin{pmatrix} I & 0 \\\ L_{1,0} & I \end{pmatrix} \begin{pmatrix} U_{0,0} & U_{0,1} \\\ 0 & S_{1,1} \end{pmatrix} $$ {{% /alert %}}
Le lien entre factorisation partielle et factorisation complète est donné par le théorème suivant :
{{< thm/thm theorem >}}
La matrice
Ce théorème nous dit que dès lors qu'on arrive à décomposer un bloc de la diagonale
Pour obtenir la factorisation complète, un algorithme itératif possible consiste à appliquer la factorisation partiellement successivement sur les compléments de Schur
{{< div class="course_lu carousel_lu" style="text-align:center">}} {{< div class="course_lu carousel-cell_lu">}} {{< svg file="course/4m053/lu_algo_0.svg" >}} {{< /div >}} {{< div class="course_lu carousel-cell_lu">}} {{< svg file="course/4m053/lu_algo_1.svg" >}} {{< /div >}} {{< div class="course_lu carousel-cell_lu">}} {{< svg file="course/4m053/lu_algo_2.svg" >}} {{< /div >}} {{< div class="course_lu carousel-cell_lu">}} {{< svg file="course/4m053/lu_algo_3.svg" >}} {{< /div >}} {{< div class="course_lu carousel-cell_lu">}} {{< svg file="course/4m053/lu_algo_4.svg" >}} {{< /div >}} {{< div class="course_lu carousel-cell_lu">}} {{< svg file="course/4m053/lu_algo_5.svg" >}} {{< /div >}} {{< /div >}}
L = 0;
U = 0;
S = A;
for k =0:N-1
// Pivot
pivot = S(k,k)
// Colonne de L
L(k,k) = 1;
for i = k+1:N-1
L(i,k) = S(i,k) / pivot;
// Ligne de U
U(k,k) = S(k,k);
for j = k+1:N-1
U(k,j) = S(k,j);
// Complément de Schur
for i = k+1:N-1
for j = k+1:N-1
S(i,j) = S(i,j) - L(i,k)*U(k,j);
Plutôt que de stocker 3 matrices L
, U
et S
, dont [on sait que cela coûte très cher]({{<relref "dense_cost.md">}}), on remarque que l'on peut se passer de ...:
- ... la matrice
S
en modifiant directementU
: le bloc$U_{k,k}$ (en "bas à droite") contiendra le complément de Schur - ... la matrice
L
en la stockant dansU
et en supprimant son terme diagonal (qui vaut 1 et peut donc devenir "implicite") - ... la matrice
U
et travailler directement dansA
Cela donne le pseudo-code suivant :
{{% div class="course_lu carousel" %}}{{% div class="course_lu carousel-cell" style="width:50%; margin-right: 10px;"%}}
S = A;
L = 0;
U = 0;
for k =0:N-1
// Pivot
pivot = S(k,k)
// Colonne de L
L(k,k) = 1;
for i = k+1:N-1
L(i,k) = S(i,k) / pivot;
// Ligne de U
U(k,k) = S(k,k);
for j = k+1:N-1
U(k,j) = S(k,j);
// Complément de Schur
for i = k+1:N-1
for j = k+1:N-1
S(i,j) = S(i,j) - L(i,k)*U(k,j);
Origine{{% /div %}}{{% div class="course_lu carousel-cell" style="width:50%; margin-right: 10px;"%}}
L = 0;
U = A;
for k =0:N-1
// Pivot
pivot = U(k,k)
// Colonne de L
L(k,k) = 1;
for i = k+1:N-1
L(i,k) = U(i,k) / pivot;
// Ligne de U
U(k,k) = U(k,k);
for j = k+1:N-1
U(k,j) = U(k,j);
// Complément de Schur
for i = k+1:N-1
for j = k+1:N-1
U(i,j) = U(i,j) - L(i,k)*U(k,j);
Suppression de S (stockée dans U){{% /div %}}{{% div class="course_lu carousel-cell" style="width:50%; margin-right: 10px;"%}}
U = A;
for k =0:N-1
// Pivot
pivot = U(k,k)
// Colonne de L
// U(k,k) = 1;
for i = k+1:N-1
U(i,k) = U(i,k) / pivot;
// Ligne de U
U(k,k) = U(k,k);
for j = k+1:N-1
U(k,j) = U(k,j);
// Complément de Schur
for i = k+1:N-1
for j = k+1:N-1
U(i,j) = U(i,j) - U(i,k)*U(k,j);
Suppression de L (stockée dans U){{% /div %}}{{% div class="course_lu carousel-cell" style="width:50%; margin-right: 10px;"%}}
for k =0:N-1
// Pivot
pivot = A(k,k)
// Colonne de L
// A(k,k) = 1;
for i = k+1:N-1
A(i,k) = A(i,k) / pivot;
// Ligne de U
A(k,k) = A(k,k);
for j = k+1:N-1
A(k,j) = A(k,j);
// Complément de Schur
for i = k+1:N-1
for j = k+1:N-1
A(i,j) -= A(i,k)*A(k,j);
Suppression de U{{% /div %}}{{% /div %}}
{{% alert exercise %}}
Avant de coder quoi que ce soit, modifiez le pseudo code de la factorisation LU
de A
effectuée directement dans la matrice A
: nettoyez le de certaines opérations rendues inutiles !
{{% /alert %}}
{{% alert exercise %}}
Implémentez une méthode de la classe Matrice
qui factorise la Matrice
sur place :
void Matrice::decomp_LU();
{{% /alert %}}
{{% alert warning %}}
Après application de l'algorithme, la Matrice A
sera modifiée de telle sorte que sa partie triangulaire inférieure soit égale à
- Le produit matrice vecteur n'a alors plus de sens une fois cet algorithme appliqué !
- Il ne faut pas ré-appliquer la factorisation LU sur A (c'est de toute façon inutile, mais une erreur arrive si vite...)
Il peut être intéressant de rajouter un paramètre à la classe Matrice
de type booleen
(un "flag") permettant de déterminer si une matrice a été, ou non, déjà factorisée.
{{% /alert %}}
{{% alert tips %}}
Une première étape pour valider votre factorisation LU : calculer le produit
{{% alert exercise %}}
Validez votre factorisation
{{% alert exercise %}}
Résolvez numériquement le problème suivant à l'aide de la factorisation LU :
$$
A X= b,
$$
où
{{% js type="text/javascript" src="https://npmcdn.com/flickity@2/dist/flickity.pkgd.js" %}} {{% js type="text/javascript" src="../lu_fact.js" %}}