# Utilisation de la décomposition $LDL^\top$ en simulation stochastique

Alban Notin, Jean-Luc Dulong, Nicolas Gayton, Pierre Villon

► To cite this version:

Alban Notin, Jean-Luc Dulong, Nicolas Gayton, Pierre Villon. Utilisation de la décomposition $LDL^\top$ en simulation stochastique. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01408464

HAL Id: hal-01408464

https://hal.science/hal-01408464v1

Submitted on 5 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Utilisation de la décomposition LDL<sup>T</sup> en simulation stochastique

A. Notin<sup>1,2</sup>, J.L. Dulong<sup>2</sup>, N. Gayton<sup>3</sup>, P. Villon<sup>2</sup>

<sup>1</sup> CETIM

52, avenue Félix Louat, B.P.80067, 60304 SENLIS Cedex, France
alban.notin@cetim.fr

<sup>2</sup> Laboratoire Roberval

UTC, CNRS FRE No. 2833, B.P.20319, 60206 Compiègne Cedex, France
{jean-luc.dulong, pierre.villon}@utc.fr

<sup>3</sup> LaMI

Université Blaise Pascal et IFMA, F-61375 Aubière, France
nicolas.gayton@ifma.fr

**Résumé** – Dans le cadre des simulations de Monte Carlo appliquées à la mécanique linéaire élastique, une alternative basée sur la décomposition LDL<sup>T</sup> de la matrice de rigidité est présentée. Cette décomposition permet de diminuer les temps de calcul des problèmes à large bande et d'accéder directement aux moments statistiques de la réponse (ici le déplacement). Ces résultats sont utilisés pour une analyse en fatigue au travers du critère de Dang Van.

**Mots clés** – calcul stochastique, décomposition LDL<sup>T</sup>, réduction de modèle, fatigue

## 1 Introduction

Ce papier introduit une utilisation originale de la décomposition LDL<sup>T</sup> dans le but d'accélérer les simulations de Monte Carlo. En effet, même si les simulations de Monte Carlo présentent de nombreux avantages [1,2], elles deviennent très exigeantes en coût de calcul et espace de stockage lorsque la taille et la complexité des problèmes augmentent. Une solution pour diminuer les temps de calcul et l'espace de stockage est d'utiliser les techniques de réduction de modèle. Ces techniques peuvent être classées en deux catégories; la première s'intéresse au modèle stochastique (utilisation par exemple du chaos polynomial [3]) et la seconde catégorie se concentre sur le modèle mécanique et sa résolution. L'approche décrite ici appartient à la deuxième catégorie de méthodes puisque qu'une nouvelle stratégie de résolution est proposée. Cette stratégie repose sur l'externalisation d'un maximum d'étapes de la boucle de calcul et a été développée dans le cadre de la mécanique linéaire élastique. Les résultats sont présentés sur une analyse à la fatigue utilisant le critère de Dang Van.

## 2 Exposé de la méthode LDL<sup>T</sup>

### 2.1 Principes généraux

Supposons une structure de module d'Young $E(x, \theta)$ modélisé par un champ stochastique et dont le chargement $f$ est déterministe. Afin de faciliter la lecture du document, le champ
stochastique $E(x, \theta)$ sera noté $E(\theta)$ dans la suite. En utilisant la formulation éléments finis en déplacement $u$ [4] dans le cadre de l'élasticité linéaire, le système linéaire à résoudre est de la forme :

$$ K(E(\theta))u(\theta) = f \qquad (1) $$

Par construction, la matrice de rigidité $K$ est symétrique, bande et définie positive. Elle admet donc une décomposition LDL$^T$ [5] :

$$ K(E(\theta)) = L(E(\theta))D(E(\theta))L^T(E(\theta)) \qquad (2) $$

$L(E(\theta))$ est une matrice triangulaire inférieure pleine, de diagonale unitaire et $D(E(\theta))$ est une matrice diagonale (dont les termes diagonaux sont notés $d_i(E(\theta))$). A partir de cette décomposition, on formule les hypothèses de travail suivantes en supposant de petites variations de $E$ :

1. les variations de la matrice $L(E(\theta))$ sont négligeables: $L(E(\theta)) = L(E(\theta_0)) = L_0$ ;

2. l'aléa est uniquement porté par la matrice diagonale $D(E(\theta))$.

$L_0$ est calculée pour une réalisation du champ stochastique notée $\theta_0 : E(\theta_0)$ est un champ déterministe fonction de l'espace et :

$$ K(E(\theta_0)) = L_0 D(E(\theta_0)) L_0^T \qquad (3) $$

Pour chaque instanciation $\theta$ du champ de module d'Young, l'objectif est donc de calculer la « meilleure matrice $D(E(\theta))$ » permettant d'approcher la matrice $K(E(\theta))$. Cette recherche de la meilleure matrice $D(E(\theta))$ se ramène donc à un problème d'optimisation qui doit être résolu pour chaque instanciation: il s'agit de minimiser l'écart entre $K(E(\theta))$ et sa décomposition LDL$^T$: $L_0 D(E(\theta)) L_0^T$. On forme pour cela la fonctionnelle $J_E(d)$ :

$$ J_E(d) = \frac{1}{2} \| K(E(\theta)) - L_0 D(E(\theta)) L_0^T \|_F^2 \qquad (4) $$

avec $\| \cdot \|_F^2$ la norme de Frobenius ($\|A\|_F^2 = \text{trace}(AA^T)$). La fonctionnelle peut donc se réécrire :

$$ J_E(d) = \frac{1}{2} \text{trace} \left\{ \left( K(E(\theta)) - L_0 D(E(\theta)) L_0^T \right) \left( K(E(\theta)) - L_0 D(E(\theta)) L_0^T \right)^T \right\} \qquad (5) $$

**Remarque 1:** La matrice $L_0 D(E) L_0^T$ peut se réécrire comme combinaison des termes $d_i$ de la matrice diagonale $D$ :

$$ L_0 D(E) L_0^T = \sum_{i=1}^{n_{ddl}} d_i(E) \Lambda^i \qquad (6) $$

Avec $n_{ddl}$ taille de la matrice et $\Lambda^i = (L_0)_{.,i} (L_0^T)_{.,i}$ (produit de la ième colonne de $L_0$ avec la ième
ligne de $L_0^T$).

**Remarque 2:** Les matrices $\Lambda^i = (L_0)_{,i} (L_0^T)_{,i}$ sont indépendantes de $E$. On les calcule une fois pour toute en dehors de la boucle de calcul de Monte Carlo.

## 2.2 Résolution du problème d'optimisation

Résoudre le problème d'optimisation (4) revient à rechercher les zéros de la dérivée de l'équation (5) pour tous les termes diagonaux $d_i$ :

$$ \left( \frac{d J_E}{d d_i} \right) \delta d_i = \operatorname{trace} \left\{ \left( K(E) - L_0 D(E) L_0^T \right) \left[ -L_0 \left[ \frac{d D(E)}{d d_i} \delta d_i \right] L_0^T \right] \right\} = 0 \quad (7) $$

ce qui équivaut à résoudre le système suivant:

$$ Ad(E) = c(E) \quad (8) $$

Avec :

$$ A_{ij} = \operatorname{trace}\{\Lambda^i \Lambda^j\} = \left[ (L_0)^T (L_0)_{,i} \right]^2 \quad (9) $$

$d(E)$ le vecteur des termes diagonaux de la matrice $D(E)$ et :

$$ c_i(E) = \operatorname{trace}\{K(E)\Lambda^i\} = (L_0^T)_{,i} K(E) (L_0)_{,i} \quad (10) $$

La matrice $A$ est une matrice bande, de même taille et largeur de bande que $K$ (la décomposition $LDL^T$ préserve la largeur de bande). Elle est indépendante des réalisations du champ $E$. De ce fait, son calcul n'est effectué qu'une seule fois en dehors de la boucle de Monte Carlo (Cf. Remarque 2).

## 2.3 Calcul des déplacements

En supposant un chargement déterministe et en injectant l'équation (3) dans l'équation (1) le déplacement s'écrit:

$$ u(E) = L_0^{-T} D^{-1}(E) L_0^{-1} f = F \left\{ \frac{1}{d(E)} \right\} \quad (11) $$

Le coût de ce calcul se résume au coût de deux remontées des matrices triangulaires bande $L_0$. La matrice $F$ est calculée hors boucle de Monte Carlo.

## 2.4 Calcul des moments statistiques

A partir de l'équation (11), les expressions de la moyenne $\bar{u}$ et de la matrice de covariance $C_u$ peuvent être directement déduites. $\bar{u}$ s'écrit:

$$ \bar{u} = L_0^{-T} \left\{ \int_{\Omega} D^{-1}(E(\theta)) dp \theta \right\} L_0^{-1} f \quad (12) $$

Et :

$$ C_u = L_0^{-T} C_d L_0^{-1} - \bar{u} \bar{u}^T \quad (13) $$
Avec :

$$ (C_d)_{ij} = \int_{\Omega} \frac{dp\theta}{d_i d_j} (C_g)_{ij} ; \quad C_g = L_0^{-1} ff^T L_0^{-T} \quad (14) $$

$\Omega$ représente l'ensemble des instanciations de $\theta$.

## 3 Application : Plaque 2D en traction pure

La méthode développée précédemment est appliquée sur une plaque en traction pure dont le module d'Young est modélisé par un champ aléatoire corrélaté 2D de moyenne $\lambda = 210000$ MPa, de coefficient de variation $\zeta$ et de modèle de corrélation [6]:

$$ R(\mathbf{x}, \mathbf{x}') = \exp\left(-\frac{|\mathbf{x}_1 - \mathbf{x}'_1|}{h_1} - \frac{|\mathbf{x}_2 - \mathbf{x}'_2|}{h_2}\right) \quad (15) $$

avec $h_1 = h_2 = 1$. La largeur $x_2$ de la plaque est de 1 et sa longueur $x_1$ vaut 2.

<figure><img src="image_2.png" /><figcaption>Figure 1 – Plaque : maillage et conditions aux limites.</figcaption></figure>

La plaque est maillée avec des éléments Q4 (figure 1). Le maillage est régulier : en modifiant le nombre d'éléments sur $x_1$ et $x_2$ il est alors possible de jouer sur le nombre total de ddl ainsi que sur la largeur de bande $b$ de la matrice de rigidité. Dans un premier temps, les calculs sont effectués sur différents maillages pour 2000 tirages aléatoires de module d'Young. Les résultats sont exprimés en temps de calcul et en FLOPS (FLoating Operations Per Second). Pour chaque tirage de Monte Carlo, l'équation (1) est résolue de manière classique par décomposition de Cholesky : pour une matrice bande, son coût est homogène à $O(n_{ddl}b^2)$. Pour la méthode LDL$^T$, le coût de résolution est dominé par le calcul de (10) qui vaut $O(n_{ddl}^2)$. La décomposition LDL$^T$ est donc moins coûteuse qu'une résolution classique si $\sqrt{n_{ddl}} < b$. Ce résultat est en accord avec les résultats de la figure 2 qui représentent l'évolution du temps de calcul pour différentes largeurs de bande en fonction de la finesse du maillage. La méthode proposée est très rapide si la structure comporte un grand nombre de ddl et a une grande largeur de bande : pour une demi-largeur de bande de 25, la simulation de Monte Carlo classique est plus performante alors que pour une largeur de 105 les deux techniques sont équivalentes. Pour une largeur de bande de 305, le coût de calcul est réduit d'un facteur 8.

Le calcul de la durée de vie s'effectue grâce au critère de Dang Van [7]. On s'intéresse à l'erreur sur les moments statistiques de la variable aléatoire « critère de Dang Van ».
L'ensemble des résultats pour ces erreurs est résumé dans le tableau 1 et les graphes de convergence pour $\zeta = 0.1$ sont présentés figure 3. L'ordre de grandeur de l'erreur est le même pour les deux coefficients de variation étudiés ($\zeta = 0.1$ et $\zeta = 0.01$). Les résultats pour la moyenne sont excellents alors que ceux pour l'écart type sont mauvais. Actuellement, nous travaillons sur l'interprétation de cet écart afin de le corriger. Une première interprétation est de considérer que le report de l'aléa uniquement sur les termes diagonaux « rigidifie » la structure. Ceci a pour conséquence de diminuer l'écart type des déplacements. Néanmoins, les résultats sur les moments d'ordre 3 et 4 suggèrent que cette rigidification n'a pas d'influence sur la forme de la densité de probabilité : on « rigidifie » sans déformer.

<figure><img src="image_3.png" /><figcaption>Figure 2 – Plaque : Comparaison des temps de calcul entre Monte Carlo et $LDL^T$ pour $b=\{25, 105, 305, 1005\}$.</figcaption></figure>

<table><caption>Table 1 - Erreur sur les moments statistiques en % en fonction de $\zeta$.</caption><thead><tr><td>Coeff. variation</td><td>Moyenne</td><td>Ecart type</td><td>Asymétrie</td><td>Aplatissement</td></tr></thead><tbody><tr><td>0.1</td><td>0.13</td><td>56.67</td><td>14.65</td><td>0.39</td></tr><tr><td>0.01</td><td>0.0011</td><td>56.02</td><td>13.31</td><td>0.24</td></tr></tbody></table>

## 4 Conclusion

L'approximation de $K$ par sa décomposition $LDL^T$ permet de s'affranchir de son inversion quand on résout (1). La matrice triangulaire inférieure $L$ est calculée une seule fois et l'inversion de $K$ se résume donc à l'inversion de $D$. Le gain en temps de calcul est proportionnel au ratio $b/\sqrt{n_{ddl}}$ et devient donc favorable pour les problèmes à large bande. Les résultats montrent le bon comportement de la méthode en termes de temps de calcul et de convergence puisque que la méthode converge à chaque fois. Des travaux complémentaires sont en cours pour corriger l'effet de « rigidification » mis en évidence.
<figure><img src="image_4.png" /><figcaption>Figure 3 – Plaque : Convergence des moments statistiques sur le critère de Dang Van pour $\zeta = 0.1$ (en pointillés, la valeur supposée exacte obtenue par simulation de Monte Carlo).</figcaption></figure>

Références

[1] M. Papadrakakis, V. Papadapoulos, Robust and efficient methods for stochastic finite element analysis using Monte Carlo simulation. *Computer Methods in Applied Mechanics and Engineering*, 134, 325–340, 1996.

[2] G. Stefanou, M. Papadrakakis, Stochastic finite element analysis of shells with combined random material and geometric properties. *Computer Methods in Applied Mechanics and Engineering*, 19, 139–160, 2004.

[3] R.G. Ghanem, P.D. Spanos, *Stochastic finite elements – A spectral approach*, Springer Verlag, 1991 (Revised Edition published by Dover in 2003).

[4] O.C. Zienkiewicz and R.L. Taylor. *The finite element method*. Butterworth-Heinemann, 2005. 6th edition

[5] P. Okunev and C. Johnson. *Necessary and sufficient conditions for existence of the LU factorization of an arbitrary matrix*. 1997.

[6] M. Lemaire. *Fiabilité des Structures - Couplage Mécano-Fiabiliste Statique*. Hermès, 2005. En collaboration avec A. Chateauneuf et J.C. Mitteau.

[7] K. Dang Van, B. Griveau, O. Message, *On a new multiaxial fatigue criterion: theory and application*. *Biaxial and Multiaxial Fatigue*. EGF 3 (Edited by M.W. Brown and K.J.Miller), London: Mechanical Engineering Publications, 1989, pp. 479-496.