# Méthode de décomposition spectrale généralisée pour la résolution des problèmes éléments finis stochastiques

Anthony Nouy

► To cite this version:

Anthony Nouy. Méthode de décomposition spectrale généralisée pour la résolution des problèmes éléments finis stochastiques. 8e Colloque national en calcul des structures, CSMA, May 2007, Giens, France. hal-01504123

HAL Id: hal-01504123

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

Submitted on 8 Apr 2017

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 CC0 1.0 - Universal - International License
# Méthode de décomposition spectrale généralisée pour la résolution des problèmes élémentsofins stochastiques

Anthony Nouy

Institut de recherche en Génie civil et mécanique (GeM),
Université de Nantes Atlantique, Ecole Centrale de Nantes, UMR CNRS 6183
anthony.nouy@univ-nantes.fr

RÉSUMÉ. Les méthodes élémentsofins stochastiques de type Galerkin conduisent généralement à la résolution de problèmes de très grande taille. Les coûts de calcul et de stockage mémoire engendrés limitent encore leur utilisation à une faible dimension stochastique. Nous proposons ici une technique d'approximation pour la résolution d'équations aux dérivées partielles stochastiques qui tente de répondre à ces difficultés. Cette technique, basée sur une décomposition spectrale généralisée de la solution, nécessite la résolution d'un faible nombre de problèmes déterministes et d'équations stochastiques. Elle conduit à la génération automatique d'une base réduite déterministe (ou stochastique), qui est quasi-optimale en regard de l'opérateur et du second membre du problème stochastique initial.

ABSTRACT. Galerkin stochastic finite element methods, when dealing with high stochastic dimension, lead to the resolution of huge systems of equations. Induced calculation costs and memory requirements still limit their use to a low stochastic dimension. Here, we propose a method which tries to overcome these difficulties. This method, based on a generalization of classical spectral decomposition techniques, needs for the resolution of a small set of uncoupled deterministic problems and stochastic equations. It leads to the automatic construction of a reduced deterministic (or stochastic) basis, which is quasi-optimal with respect to the operator and right-hand side of the stochastic problem.

MOTS-CLÉS : Méthodes de calcul stochastique, Eléments finis stochastiques, Karhunen-Loeve, Décomposition spectrale généralisée

KEYWORDS: Computational stochastic methods, Stochastic finite elements, Karhunen-Loeve, Generalized Spectral Decomposition
# 1. Introduction

Les méthodes éléments finis stochastiques récemment proposées pour la résolution des équations aux dérivées partielles stochastiques constituent un outil essentiel pour la prise en compte de l'aléa dans le calcul de structure. Les méthodes de type Galerkin (Ghanem et al., 1991, Deb et al., 2001, Matthies et al., 2005) conduisent de manière systématique à une solution d'une grande précision, explicite en fonction des variables aléatoires décrivant l'aléa. Cependant, le coût de calcul engendré et les ressources mémoire requises limitent encore leur utilisation à une faible dimension stochastique lorsqu'on considère des problèmes de grande taille au niveau déterministe.

Nous présentons ici une méthode alternative pour la résolution de problèmes stochastiques, inspirée d'une méthode de résolution d'équations d'évolution (Nouy et al., 2004). Cette méthode, qui généralise la technique de décomposition spectrale classique, conserve les avantages des approches de Galerkin tout en diminuant drastiquement les coûts de calcul et de stockage mémoire.

# 2. Formulation du problème stochastique et approximation

On considère que la modélisation probabiliste du problème a conduit à la définition d'un espace probabilisé adapté $(\Theta, \mathcal{B}, P)$ de dimension finie, où $\theta \in \Theta$ désigne un événement élémentaire, $\mathcal{B}$ une $\sigma$-algèbre sur $\Theta$ et $P$ une mesure de probabilité. On définit alors un espace ad-hoc $S$ de variables aléatoires définies sur $\Theta$ et à valeur dans $\mathbb{R}$. Les méthodes de Galerkin au niveau stochastique consistent à choisir un espace d'approximation $S_P = \{v = \sum_{\alpha \in \mathcal{I}_P} v_\alpha H_\alpha(\theta), v_\alpha \in \mathbb{R}, H_\alpha \in S\}$ de dimension finie. On utilise classiquement une approximation par chaos polynômial (Ghanem et al., 1991, Soize et al., 2004), où la base d'approximation $\{H_\alpha\}$ est constituée de polynômes orthogonaux pour le produit scalaire naturel sur $S$, ou encore une approximation polynômiale par morceaux (Deb et al., 2001, Le Maître et al., 2004). Un problème stochastique discrétisé au niveau déterministe et stochastique s'écrira alors généralement sous la forme : trouver $\mathbf{u} \in \mathbb{R}^n \otimes S_P$ telle que

$$E(\mathbf{v}^T \mathbf{A}(\mathbf{u})) = E(\mathbf{v}^T \mathbf{b}) \quad \forall \mathbf{v} \in \mathbb{R}^n \otimes S_P \quad [1]$$

où $\mathbf{b}: \Theta \to \mathbb{R}^n$ est un vecteur aléatoire et, dans le cas de problèmes linéaires, $\mathbf{A}: \Theta \to \mathbb{R}^{n \times n}$ est une matrice aléatoire. Le système [1] est un système de très grande taille $n \times P$, où $P$ est le nombre de fonctions de base stochastiques. Les techniques de résolution classiques sont de type Krylov (gradient conjugué préconditionné, conjugate gradient square, ...) (Pellissetti et al., 2000, Keese et al., 2005). La difficulté pour construire des préconditionneurs efficaces et les ressources mémoire requises par ces techniques limitent encore leur utilisation à une faible dimension stochastique lorsqu'on considère des problèmes de grande taille au niveau déterministe.
### 3. Décomposition spectrale généralisée

La démarche proposée ici consiste à approcher la solution $\mathbf{u}$ du problème [1] par

$$ \mathbf{u}^{(M)}(\theta) = \sum_{i=1}^{M} \lambda_i(\theta) \mathbf{U}_i \quad [2] $$

où les $\lambda_i \in S_P$ sont des variables aléatoires (ou fonctions stochastiques) et les $\mathbf{U}_i \in \mathbb{R}^n$ des vecteurs déterministes. Une décomposition de ce type est dite optimale si le nombre de termes $M$ est minimal pour une qualité donnée de l'approximation. L'ensemble des fonctions déterministes (resp. stochastiques) peuvent alors être considérées comme une base réduite déterministe (resp. stochastique) optimale. Ici, ni les $\lambda_i$ ni les $\mathbf{U}_i$ ne sont fixées a priori. Les questions clés sont donc : comment définir la base réduite "optimale" et comment la calculer ? En fait, la décomposition obtenue dépend de ce qu'on entend par "optimal". Si la solution était connue, la meilleure approximation serait classiquement définie en minimisant la distance à la solution au sens d'une norme sur $\mathbb{R}^n \otimes S$. Il est bien connu que la meilleure approximation $\mathbf{u}^{(M)}$ au sens de la norme naturelle sur $L^2(\Theta, dP; \mathbb{R}^n)$ est la décomposition de Karhunen-Loeve (décomposition spectrale classique) tronquée à l'ordre $M$ (Matthies et al., 2005).

Ici, nous proposons une manière simple et intuitive de définir la meilleure décomposition qui généralise la technique de décomposition spectrale classique. En notant $\mathbf{W} = (\mathbf{U}_1 \dots \mathbf{U}_M) \in \mathbb{R}^{n \times M}$ et $\mathbf{\Lambda} = (\lambda_1 \dots \lambda_M)^T \in \mathbb{R}^M \otimes S_P$, l'approximation recherchée s'écrit $\mathbf{u} = \mathbf{W}\mathbf{\Lambda}$. On la définit par les deux propriétés suivantes :

$$ E(\mathbf{\Lambda}^T (\widetilde{\mathbf{W}}^T \mathbf{A}(\mathbf{W}\mathbf{\Lambda}))) = E(\mathbf{\Lambda}^T \widetilde{\mathbf{W}}^T \mathbf{b}) \quad \forall \widetilde{\mathbf{W}} \in \mathbb{R}^{n \times M} \quad [3] $$

$$ E(\tilde{\mathbf{\Lambda}}^T (\mathbf{W}^T \mathbf{A}(\mathbf{W}\mathbf{\Lambda}))) = E(\tilde{\mathbf{\Lambda}}^T \mathbf{W}^T \mathbf{b}) \quad \forall \tilde{\mathbf{\Lambda}} \in \mathbb{R}^M \otimes S_P \quad [4] $$

Si $\mathbf{\Lambda}$ (resp. $\mathbf{W}$) était fixé, l'équation [3] (resp. [4]) définirait l'approximation de Galerkin naturelle. Ici, aucune des fonctions n'étant fixée a priori, on cherchera naturellement des fonctions vérifiant ces équations simultanément. Le système d'équations ([3],[4]) n'a pas une solution unique. Dans le cas de problèmes symétriques, s'écrivant comme un problème de minimisation d'une fonctionnelle, le meilleur choix de fonctions $(\mathbf{W}, \mathbf{\Lambda})$ est naturellement celui minimisant la fonctionnelle. Dans ce cas, les équations [3] et [4] expriment simplement les conditions de stationnarité de la fonctionnelle. On montre que cette définition conduit à la résolution d'un problème aux valeurs propres "étendu", la meilleure approximation étant celle associée aux $M$ plus grandes valeurs propres "généralisées" (Nouy, 2007). Pour construire cette approximation, il reste donc à construire des algorithmes capables de capter efficacement le haut de ce spectre "généralisé".

### 4. Algorithme de construction de l'approximation

Un algorithme de type méthode des puissances a déjà été proposé (Nouy et al., 2006, Nouy, 2007). Il consiste à construire successivement les couples $(\mathbf{U}_i, \lambda_i)$,
chaque couple étant obtenu en résolvant de manière itérative le système d'équations [3] et [4]. L'obtention de chaque couple ne nécessite donc que la résolution de systèmes d'équations déterministes de taille *n* et d'équations stochastiques (systèmes de taille *P*). Dans le cas d'un opérateur **A** déterministe, on montre que le problème aux valeurs propres "étendu" se dégénère en un problème aux valeurs propres classique, où la meilleure approximation correspond à une décomposition spectrale de **b** dans la métrique induite par **A**. L'algorithme proposé correspond alors à une méthode des puissances classique permettant bien de capter les plus grandes valeurs propres. Cet algorithme, bien que conduisant à une approximation satisfaisante, peut cependant s'avérer moins efficace (convergence lente de la méthode des puissances) dans le cas général, notamment pour des problèmes non-symétriques. On propose ici un nouvel algorithme (A-GSD), inspiré de la technique d'Arnoldi pour la résolution de problèmes aux valeurs propres (Saad, 1992).

**Algorithm 1 Arnoldy-type Generalized Spectral Decomposition (A-GSD)**

1: Initialize $\lambda \in S_P$

2: **for** $i = 1$ to $M_{max}$ **do**

3: Solve a deterministic problem : $\tilde{U}_i := E(\mathbf{A}\lambda^2)^{-1}E(\mathbf{b}\lambda)$

4: Orthogonalize : $\mathbf{U}_i := \tilde{U}_i - \sum_{j=1}^{i-1} \mathbf{U}_j^T \tilde{U}_i$

5: **if** $||\mathbf{U}_i|| < \epsilon||\tilde{\mathbf{U}}_i||$ **then break**
**else** Normalize $\mathbf{U}_i$, set $\mathbf{W} = (\mathbf{U}_1 \dots \mathbf{U}_i)$ **endif**

6: Update $\lambda$ by solving : $E(\tilde{\lambda}(\mathbf{U}_i^T \mathbf{A} \mathbf{U}_i)\lambda) = E(\tilde{\lambda}(\mathbf{U}_i^T \mathbf{b}))$, $\forall \tilde{\lambda} \in S_P$

7: **end for**

8: Find $\mathbf{\Lambda}$ by solving equation [4]

9: Update solution by adding $\mathbf{W}\mathbf{\Lambda}$, Update residual $\mathbf{b}$ and restart at step 1.

**5. Exemple : problème de thermique instationnaire**

On considère un problème de thermique instationnaire défini sur le domaine $\Omega$ (figure 1), et l'intervalle de temps $(0, 3.10^4)$. Des conditions homogènes de Dirichlet sont imposées sur $\partial\Omega_1$. Un flux de chaleur normal $g_1$ (resp. $g_2$) est imposé sur $\Gamma_1$ (resp. $\Gamma_2$). On prendra $g_1(s, \theta) = \xi_1(\theta) = \mathcal{N}(10^4, 0.2 \ 10^4)$ et $g_2(s, \theta) = \xi_2(\theta) = \mathcal{N}(10^4, 0.2 \ 10^4)$, des variables aléatoires gaussiennes<fnref n="1" />. Le problème est formulé sur $u = T - T_0$ où $T$ est la température et $T_0$ la température initiale supposée homogène. Il se met sous la forme variationnelle suivante :

$$u \in \mathcal{V}_n \otimes S_P, \quad E(a(u, v; \theta)) = E(b(v; \theta)) \quad \forall v \in \mathcal{V}_n \otimes S_P$$

$$a(u, v) = \int_{0}^{T} \int_{\Omega} \left( c \frac{\partial u}{\partial t} v + \kappa \nabla u \cdot \nabla v \right) dx dt, \quad b(v) = \sum_{i=1}^{2} \int_{0}^{T} \int_{\Gamma_i} g_i v \, ds \, dt$$

où $c$ et $\kappa$ sont les paramètres matériau, pris comme deux variables aléatoires uniformes<fnref n="2" /> : $c = \xi_3(\theta) = U(2.34 \ 10^6, 4.68 \ 10^6)$, $\kappa = \xi_4(\theta) = U(30, 60)$. Les variables

<footnote n="1">1. $\mathcal{N}(\mu, \sigma)$ désigne une variable gaussienne de moyenne $\mu$ et d'écart-type $\sigma$.</footnote>

<footnote n="2">2. $U(\alpha, \beta)$ désigne une variable uniforme sur $[\alpha, \beta]$.</footnote>
aléatoires {$\xi_i$}$_{i=1}^4$ sont considérées indépendantes. On choisit au niveau stochastique une approximation par chaos polynômial généralisé de degré $p$ en dimension $m=4$. L'espace d'approximation $S_P$ est donc de dimension $P = \frac{(m+p)!}{m!p!}$. Au niveau spatial, on utilise une approximation élément fini ($n_s = 1179$ degrés de liberté) et au niveau temporel, une approximation de Galerkin discontinue ($n_t = 50$ intervalles temporels et approximation constante par intervalle). L'espace d'approximation $V_n$ est de dimension $n = n_s \times n_t$. On notera que le problème discrétisé peut encore s'écrire sous la forme [1]. Cependant, les équations déterministes de type [3] sont bien sûr résolues de manière incrémentale en temps (intervalle par intervalle). Une solution de référence est calculée par une résolution incrémentale en temps associée à un gradient conjugué préconditionné (PCG) pour le problème stochastique défini sur chaque intervalle temporel. La figure 2 montre l'erreur de l'approximation $\mathbf{u}^{(M)}$ calculée par l'algorithme (A-GSD) (norme sur $L^2(0, T; L^2(\Omega)) \otimes L^2(\Theta)$). On observe une bonne convergence en fonction de l'ordre de la décomposition $M$ (seulement 4 fonctions suffisent à obtenir une erreur globale de $10^{-2}$). L'approche semble relativement robuste. On notera notamment que dans une approche classique, il est difficile d'obtenir une telle erreur globale en raison du cumul d'erreur au niveau temporel. Pour une erreur de $10^{-2}$ (M=4), un facteur 176 est gagné sur le stockage mémoire de la solution. Le tableau 1 compare les temps de calcul en fonction de la discrétisation déterministe et stochastique. On observe qu'à $n$ fixé, la dimension stochastique a très peu d'influence, contrairement à une approche classique (PCG). Sur ce problème très simple, on peut déjà gagner jusqu'à un facteur 22 sur le temps de calcul.

<figure><img src="image_2.png" /><figcaption>Figure 1. L-shape domain</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2. Convergence de la décomposition $\mathbf{u}^{(M)}$</figcaption></figure>

## 6. Conclusion

On a proposé une méthode robuste et générale pour la résolution d'équations stochastiques. Cette méthode permet la construction *a priori* d'une décomposition spectrale quasi-optimale de la solution, et conduit à des gains considérables en temps de calcul et stockage mémoire. Les travaux futurs concernent l'extension de cette méthode aux problèmes non-linéaires.
<table><caption>**Tableau 1.** *Comparaison A-GSD / PCG : temps de calcul (s) en fonction de n et P*</caption><thead><tr><th rowspan="2">P</th><th colspan="2">n = 690<br>(69 × 10)</th><th colspan="2">n = 8670<br>(289 × 30)</th><th colspan="2">n=26360<br>(659 × 40)</th><th colspan="2">n=58950<br>(1179 × 50)</th></tr><tr><th>A-GSD</th><th>PCG</th><th>A-GSD</th><th>PCG</th><th>A-GSD</th><th>PCG</th><th>A-GSD</th><th>PCG</th></tr></thead><tbody><tr><td>35</td><td>0.7</td><td>0.47</td><td>2.2</td><td>6</td><td>10.2</td><td>18</td><td>33</td><td>49</td></tr><tr><td>126</td><td>0.7</td><td>1.1</td><td>2.3</td><td>16</td><td>10.2</td><td>66</td><td>33</td><td>150</td></tr><tr><td>330</td><td>1</td><td>3.1</td><td>2.8</td><td>48</td><td>10.4</td><td>160</td><td>33</td><td>355</td></tr><tr><td>715</td><td>2.2</td><td>13.4</td><td>3.8</td><td>106</td><td>11.3</td><td>323</td><td>34</td><td>742</td></tr></tbody></table>

Remerciements

Ces travaux sont soutenus par l'Agence National de la Recherche (ANR) dans le cadre du programme "Jeunes chercheurs" (projet "X-SFEM").

**7. Bibliographie**

Deb M., Babuška I., Oden J. T., « Solution of stochastic partial differential equations using Galerkin finite element techniques », *Comp. Meth. App. Mech. Eng.*, vol. 190, p. 6359-6372, 2001.

Ghanem R., Spanos P., *Stochastic finite elements : a spectral approach*, Springer, Berlin, 1991.

Keese A., Mathhies H. G., « Hierarchical parallelisation for the solution of stochastic finite element equations », *Comp. Meth. App. Mech. Eng.*, vol. 83, p. 1033-1047, 2005.

Le Maître O. P., Knio O. M., Najm H. N., Ghanem R. G., « Uncertainty propagation using Wiener-Haar expansions », *Journal of Computational Physics*, vol. 197, p. 28-57, 2004.

Matthies H. G., Keese A., « Galerkin methods for linear and nonlinear elliptic stochastic partial differential equations », *Comp. Meth. App. Mech. Eng.*, vol. 194, n° 12-16, p. 1295-1331, 2005.

Nouy A., « A generalized spectral decomposition technique to solve stochastic partial differential equations », *Submitted to Comp. Meth. App. Mech. Eng.*, 2007.

Nouy A., Ladevèze P., « Multiscale computational strategy with time and space homogenization : a radial-type approximation technique for solving micro problems », *International Journal for Multiscale Computational Engineering*, vol. 170, n° 2, p. 557-574, 2004.

Nouy A., Schoefs F., « A generalized spectral decomposition technique to solve stochastic partial differential equations », *7th World Congress on Computational Mechanics (WCCM VII), Los Angeles, USA*, 2006.

Pellissetti M. F., Ghanem R. G., « Iterative solution of systems of linear equations arising in the context of stochastic finite elements », *Adv. Eng. Soft.*, vol. 31, p. 607-616, 2000.

Saad Y., *Numerical methods for large eigenvalue problems*, Halstead Press, New York, 1992.

Soize C., Ghanem R., « Physical systems with random uncertainties : chaos representations with arbitrary probability measure », *SIAM J. Sci. Comput.*, vol. 26, n° 2, p. 395-410, 2004.