3. Principe de la résolution par élimination#

3.1. Définition du problème à résoudre#

Le système à résoudre peut se mettre sous la forme générale suivante:

(3.1)#\[\left[K\right]\lbrace u\rbrace =\lbrace f\rbrace\]

sous la contrainte

(3.2)#\[\left[B\right]\lbrace u\rbrace =\lbrace {u}_{0}\rbrace\]

Avec l’introduction de multiplicateurs de Lagrange, le problème devient

(3.3)#\[\begin{split}\lbrace \begin{array}{c}\left[K\right]\lbrace u\rbrace +{\left[B\right]}^{T}\lbrace \lambda \rbrace =\lbrace f\rbrace \\ \left[B\right]\lbrace u\rbrace =\lbrace {u}_{0}\rbrace \end{array}\end{split}\]

et se met sous la forme générale

(3.4)#\[\begin{split}\left[\begin{array}{cc}K& {B}^{T}\\ B& 0\end{array}\right]\left\lbrace \begin{array}{c}u\\ \lambda \end{array}\right\rbrace =\left\lbrace \begin{array}{c}f\\ {u}_{0}\end{array}\right\rbrace\end{split}\]

Dans le cas de la double dualisation (AFFE_CHAR_MECA), cette matrice est augmentée, et le système à résoudre devient

(3.5)#\[\begin{split}\left[\begin{array}{ccc}K& {B}^{T}& {B}^{T}\\ B& -\alpha {I}_{d}& \alpha {I}_{d}\\ B& \alpha {I}_{d}& -\alpha {I}_{d}\end{array}\right]\left\lbrace \begin{array}{c}u\\ {\lambda}_{1}\\ {\lambda}_{2}\end{array}\right\rbrace =\left\lbrace \begin{array}{c}f\\ {u}_{0}\\ {u}_{0}\end{array}\right\rbrace\end{split}\]

après avoir posé naturellement \({\lambda}_{1}={\lambda}_{2}=\lambda /2\) .

Pour rendre son caractère semi-défini positif au problème générique \(\mathit{Ku}=f\) , dans le cas où la matrice \(K\) présente une topologie découlant de la mise sous la forme () ou (), on va chercher à décomposer la solution \(u\) du problème sous la forme \(u={u}_{p}+{u}_{g}\) , où

  • \({u}_{p}\) est la solution particulière du problème,

  • \({u}_{g}\) est la solution générale du problème, associée aux conditions aux limites \([B]\lbrace {u}_{g}\rbrace =\lbrace 0\rbrace\) .

3.2. Recherche de la solution particulière#

La recherche de la solution particulière consiste simplement à trouver \({u}_{p}\) tel que \([B]\lbrace {u}_{p}\rbrace =\lbrace {u}_{0}\rbrace\) . Pour calculer \({u}_{p}\) , on considère la décomposition sur l’image et le noyau de \(B\) , soit \({u}_{p}={u}_{p-i}+{u}_{p-n}\) , avec \({u}_{p-i}\) appartenant à l’image de \(B\) et \({u}_{p-n}\) appartenant au noyau de \(B\) .

Ce problème peut être résolu au sens des moindres carrés, et dans ce cas, on considère la solution de () de norme \({L}_{2}\) minimale, soit \({u}_{p}={u}_{p-i}\) . \(\lbrace {u}_{p-i}\rbrace\) est la solution du problème de minimisation:

(3.6)#\[\lbrace {u}_{p-i}\rbrace =\underset{\lbrace {u}_{p}\rbrace }{\text{ArgMin}}\left({\Vert [B]\lbrace {u}_{p}\rbrace -\lbrace {u}_{0}\rbrace \Vert }^{2}\right)\]

Pour calculer efficacement \({u}_{p-i}\) , sous réserve que \(B\) soit de rang plein (rang égal à \(\mathit{Nb}\) , donc pas de contraintes redondantes), on exploite la propriété suivante:

(3.7)#\[\left([B]{[B]}^{T}\right){\left([B]{[B]}^{T}\right)}^{-1}={\left[{I}_{d}\right]}_{\mathit{NbxNb}}\]

et donc

(3.8)#\[\left([B]{[B]}^{T}\right){\left([B]{[B]}^{T}\right)}^{-1}\lbrace {u}_{0}\rbrace =\lbrace {u}_{0}\rbrace\]

La solution recherchée est naturellement donnée par:

(3.9)#\[\lbrace {u}_{p}\rbrace =\lbrace {u}_{p-i}\rbrace ={[B]}^{T}{\left([B]{[B]}^{T}\right)}^{-1}\lbrace {u}_{0}\rbrace\]

En pratique, on effectue une factorisation de Cholesky de la matrice symétrique \({\mathit{BB}}^{T}\) , à l’aide de la bibliothèque MUMPS, puis on résout

(3.10)#\[[B]{[B]}^{T}\lbrace y\rbrace =\lbrace {u}_{0}\rbrace\]

Enfin, on obtient \(\lbrace {u}_{p}\rbrace\) :

\(\lbrace {u}_{p}\rbrace =\lbrace {u}_{p-i}\rbrace ={[B]}^{T}\lbrace y\rbrace\)

3.3. Recherche de la solution générale#

La solution générale \({u}_{g}\) vérifiant \([B]\lbrace {u}_{g}\rbrace =\lbrace 0\rbrace\) , elle appartient au noyau de \(B\) . La recherche de la solution générale passe donc par une étape de construction d’une base \(Z\) du noyau de \(B\) . On peut ensuite projeter le problème sur le noyau, et résoudre ce problème, mieux posé et de taille réduite.

3.3.1. Obtention d’une base du noyau#

Formellement, la construction du noyau de \(B\) peut se faire de plusieurs façons:

  • Décomposition en valeurs singulières de \(B\) , soit \(B={\mathit{USV}}^{T}\) . \(\mathit{Ker}(B)\) est l’espace engendré par les colonnes de \({V}^{T}\) associées à une valeur singulière nulle.

  • Décomposition \(\mathit{QR}\) de \({B}^{T}\) , soit \(B={R}^{T}{Q}^{T}\) . \(\mathit{Ker}(B)\) est l’espace engendré par les colonnes de \(Q\) associées à une valeur nulle sur la diagonale de \({R}^{T}\) .

  • Décomposition \(\mathit{LU}\) de \({B}^{T}\) . Le calcul du noyau n’est pas direct, et s’obtient à partir de sous-blocs de \(L\) .

C’est cette troisième méthode qui est mise en œuvre dans code_aster. On utilise la bibliothèque SuperLU pour construire la factorisation \(\mathit{LU}\) de \({B}^{T}\) :

\[ \begin{align}\begin{aligned}:label: eq-11\\ {P}_{r}{B}^{T}{P}_{c}^{T}=\left[L\right]\left[U\right]=\left[\begin{array}{cc}{[{L}_{1}]}_{{N}_{b}\times {N}_{b}}& [:ref:\\Les matrices :math:`{P}_{r}` et :math:`{P}_{c}` désignent respectivement les matrices de permutation imposées aux lignes et aux colonnes de :math:`{B}^{T}` pour la construction de la factorisation :math:`\mathit{LU}` . Les blocs :math:`{L}_{1}` et :math:`{L}_{3}` sont triangulaires inférieurs et le bloc :math:`{U}_{1}` est triangulaire supérieur.\\Remarque:\\En pratique, :math:`{P}_{c}=\mathit{Id}` , il n’y a pas de permutation des colonnes de :math:`{B}^{T}` .\\On en déduit l’expression suivante pour :math:`B` :\end{aligned}\end{align} \]
(3.11)#\[B={U}^{T}{L}^{T}{P}_{r}=\left[\begin{array}{cc}{[{U}_{1}]}_{{N}_{b}\times {N}_{b}}^{T}& {[:ref:\]

Il est ensuite facile de voir qu’une base de \(\mathit{Ker}(B)\) est donnée par:

(3.12)#\[\begin{split}Z={P}_{r}^{-1}\left[\begin{array}{c}-{[{L}_{1}]}^{-T}{[{L}_{2}]}^{T}\\ {[{I}_{d}]}_{(N-{N}_{b})\times (N-{N}_{b})}\end{array}\right]\end{split}\]

On fait ici l’hypothèse que \(B\) est de rang plein. Si tel n’est pas le cas, alors il suffit de construire \({L}_{1}\) à partir du bloc de \(L\) associé aux termes non nuls de la diagonale de \(U\) .

3.3.2. Projection et résolution#

Une fois que la base du noyau est calculée, la résolution est directe. On pose \({u}_{g}=Z\overline{u}\) , et la solution du problème est donc directement donnée par

(3.13)#\[{[Z]}^{T}[K][Z]\lbrace \overline{u}\rbrace ={Z}^{T}\lbrace f\rbrace\]

puisque, par construction,

(3.14)#\[{[Z]}^{T}{[B]}^{T}\lbrace \lambda \rbrace =\lbrace 0\rbrace\]

Dans certains cas, il est nécessaire d’accéder aux multiplicateurs de Lagrange. Il faut donc les reconstruire, à partir des solutions \({u}_{p}\) et \({u}_{g}\) . Pour ce faire, on doit résoudre la première ligne du système (), soit

(3.15)#\[{[B]}^{T}\lbrace \lambda \rbrace =\lbrace f\rbrace -[K]\left([Z]\lbrace \overline{u}\rbrace +\lbrace {u}_{p}\rbrace \right)\]

Si \(B\) est de rang plein, la solution s’écrit formellement:

(3.16)#\[\lbrace \lambda \rbrace ={\left([B]{[B]}^{T}\right)}^{-1}[B]\lbrace f\rbrace -[K]\left([T]\lbrace \overline{u}\rbrace +\lbrace {u}_{p}\rbrace \right)\]

et se calcule en pratique en réutilisant la factorisationde \(B{B}^{T}\) calculée pour la recherche de la solution particulière. De façon analogue, si \(B\) n’est pas de rang plein, il est possible de supprimer les contraintes redondantes.