# Méthode de construction de bases spectrales généralisées pour l'approximation de problèmes stochastiques

Anthony Nouy

► To cite this version:

Anthony Nouy. Méthode de construction de bases spectrales généralisées pour l'approximation de problèmes stochastiques. 18e congrès de mécanique, Aug 2007, Grenoble, France. hal-01379708

HAL Id: hal-01379708

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

Submitted on 12 Oct 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
# Méthode de construction de bases spectrales généralisées pour l'approximation de problèmes 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

## Résumé :

Les méthodes éléments finis 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) dont on peut prouver l'optimalité vis-à-vis de l'opérateur et du second membre du problème stochastique.

## Abstract :

Galerkin stochastic finite elements 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 technique, 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 can be shown to be optimal with respect to the operator and right-hand side of the stochastic problem.

## Mots-clefs :

stochastic finite elements, generalized spectral decomposition, reduced basis

## 1 Introduction

Les méthodes éléments finis stochastiques offrent aujourd'hui un outil robuste pour la résolution de problèmes stochastiques. Les approches de type Galerkin au niveau stochastique [Ghanem et al. (1991); Deb et al. (2001); Matthies et al. (2005)] permettent d'obtenir une solution d'une grande qualité. La recherche de cette solution conduit cependant à résoudre des systèmes de très grande taille. Les coûts de calcul (en temps et en stockage mémoire) engendrés par les techniques de résolution classiques de type Krylov [Ghanem et al. (1996); Pellissetti et al. (2000); Keese et al. (2005)] limitent cependant l'utilisation de ces approches à une faible dimension au niveau stochastique.

Une méthode alternative pour la résolution de problèmes éléments finis stochastiques a récemment été proposée [Nouy (2007)]. La solution y est décomposée sous la forme d'une somme de fonctions, chaque fonction étant le produit d'une fonction déterministe et d'une fonction aléatoire scalaire. Contrairement aux approches classiques, ni les fonctions déterministes ni les fonctions aléatoires ne sont fixées a priori. Elles sont définies de telle manière à minimiser le nombre de fonctions nécessaires dans la décomposition pour une précision donnée. On montre alors que les fonctions sont solutions d'un problème qui peut être interprété comme un problème aux valeurs propres "étendu". L'approximation obtenue peut ainsi être interprétée comme une décomposition spectrale généralisée de la solution. Cette méthode généralise en quelque sorte la méthode de Karhunen-Loève à la résolution d'un problème stochastique. Elle
permet une réduction importante des coûts de calcul par rapport aux approches classiques. De plus, elle conduit naturellement à la construction d'une base réduite de fonctions déterministes (ou de fonctions stochastiques) qui peut être judicieusement exploitée dans le cadre de multi-résolutions.

Dans la partie 2, nous rappelons la formulation et l'approximation classique d'un problème stochastique. Dans la partie 3, nous introduisons les bases de la méthode de décomposition spectrale généralisée ainsi qu'un algorithme de construction adapté. Nous évoquons dans la partie 4 l'intérêt d'une réutilisation des bases spectrales dans le cadre de problèmes non-linéaires. Enfin, nous illustrons la méthode dans la partie 5 en l'appliquant à un problème de thermique non-linéaire.

## 2 Formulation et approximation d'un problème stochastique

La modélisation probabiliste d'un problème de physique se ramène généralement à identifier les sources d'aléas en fonction d'un nombre fini de variables aléatoires. Cela conduit à définir un espace probabilisé de dimension finie adapté à l'étude, noté $(\Theta, \mathcal{B}, P)$, où $\Theta$ désigne l'espace des événements élémentaires, $\mathcal{B}$ une tribu d'événements et $P$ une mesure de probabilité. Un élément $\theta \in \Theta$ désigne un événement élémentaire. Une première étape d'approximation au niveau déterministe (par éléments finis, différences finies, ...) permet de formuler un problème linéaire de physique sous la forme faible semi-discrétisée suivante : trouver un vecteur aléatoire $\mathbf{u} \in \mathbb{R}^n \otimes \mathcal{S}$, où $\mathcal{S}$ désigne un espace de variables aléatoires ad hoc, vérifiant

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

où $\mathbf{b} : \Theta \to \mathbb{R}^n$ est un vecteur aléatoire et où $\mathbf{A} : \Theta \to \mathbb{R}^{n \times n}$ est une matrice aléatoire.

L'approximation au niveau stochastique consiste à définir un espace d'approximation $\mathcal{S}_P = \{v = \sum_{\alpha \in \mathcal{J}_P} v_\alpha H_\alpha(\theta), v_\alpha \in \mathbb{R}, H_\alpha \in \mathcal{S}\}$ de dimension finie $P$. 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 $\mathcal{S}$, ou encore une approximation polynômiale par morceaux [Deb et al. (2001); Le Maître et al. (2004)]. L'approximation de Galerkin $\mathbf{u} \in \mathbb{R}^n \otimes \mathcal{S}_P$ est alors déterminée par la résolution du système d'équations linéaires suivant :

$$\sum_{\alpha \in \mathcal{J}_P} E(\mathbf{A} H_\alpha H_\beta) \mathbf{u}_\alpha = E(H_\beta \mathbf{b}) \quad \forall \beta \in \mathcal{J}_P \quad (2)$$

Le système (2) est un système de très grande taille $n \times P$. Les techniques de résolution classiquement utilisées sont des techniques itératives de type Krylov (gradient conjugué préconditionné (PCG), conjugate gradient square (CGS), ...) [Ghanem et al. (1996); Pellissetti et al. (2000); Keese et al. (2005)]. Ces algorithmes profitent de la structure creuse du système. Cependant, les coûts de résolution engendrés limitent l'utilisation de ces techniques à une faible dimension au niveau stochastique. De plus, le coût de stockage mémoire de la solution (et éventuellement des sous espaces de Krylov) devient considérable dans le cas de problèmes de grande dimension.

La technique proposée par la suite permet de pallier ces inconvénients tout en conservant la robustesse et la précision des approches de Galerkin au niveau stochastique.

## 3 Méthode de décomposition spectrale généralisée

### 3.1 Principe

L'idée de la méthode de décomposition spectrale généralisée (GSD) consiste à chercher une approximation du problème (1) sous la forme

$$\mathbf{u}(\theta) \approx \sum_{i=1}^{M} \lambda_i(\theta) \mathbf{U}_i \quad (3)$$
où les $\lambda_i \in S_P$ sont des variables aléatoires et les $U_i \in \mathbb{R}^n$ des vecteurs déterministes. Une décomposition de ce type sera dite optimale si le nombre de termes $M$ est minimum pour une qualité donnée de l'approximation. L'ensemble des vecteurs déterministes (resp. des variables aléatoires) peut alors être considéré comme une base réduite optimale déterministe (resp. stochastique). Un premier intérêt de cette décomposition apparaît immédiatement : le coût de stockage mémoire d'une solution est divisé par un facteur $\frac{n \times P}{M(n+P)}$. Si $M$ reste faible, le gain devient donc considérable quand les dimensions $n$ et $P$ augmentent.

Deux questions clés se posent alors naturellement : comment définir la décomposition "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 $u^{(M)}$ au sens de la norme naturelle sur $L^2(\Theta, dP; \mathbb{R}^n)$ est la décomposition de Karhunen-Loeve (ou décomposition spectrale classique) tronquée à l'ordre $M$ [Matthies et al. (2005)]. Bien sûr, la solution $u$ n'étant pas connue a priori, ce type de décomposition n'est pas accessible directement.

Nous introduisons par la suite une manière simple et intuitive de définir a priori cette décomposition ainsi qu'un algorithme de construction adapté. La méthode qui en découle généralise la méthode de décomposition spectrale classique.

## 3.2 Définition et construction de la décomposition spectrale généralisée

En notant $\mathbf{W} = (U_1 \dots U_M) \in \mathbb{R}^{n \times M}$ et $\mathbf{\Lambda} = (\lambda_1 \dots \lambda_M)^T \in \mathbb{R}^M \otimes S_P$, la décomposition recherchée s'écrit $\mathbf{u} = \mathbf{W}\mathbf{\Lambda}$. On choisit alors de définir la "meilleure" approximation de ce type par les deux propriétés d'orthogonalité suivantes :

$$E(\mathbf{\Lambda}^T(\tilde{\mathbf{W}}^T\mathbf{A}\mathbf{W}\mathbf{\Lambda})) = E(\mathbf{\Lambda}^T\tilde{\mathbf{W}}^T\mathbf{b}) \quad \forall \tilde{\mathbf{W}} \in \mathbb{R}^{n \times M} \quad (4)$$

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

Si $\mathbf{\Lambda}$ (resp. $\mathbf{W}$) était fixé, l'équation (4) (resp. (5)) 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. 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 (4) et (5) expriment simplement les conditions de stationnarité de la fonctionnelle. On montre alors 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" (cf. [Nouy (2007)]). L'approximation obtenue est dans ce cas la meilleure approximation du type (3) au sens de l'opérateur et du second membre.

Pour construire l'approximation (3), il reste à proposer un algorithme capable de capter efficacement le haut de ce "spectre généralisé". L'algorithme 1 proposé est inspiré de la méthode des puissances pour la résolution des problèmes aux valeurs propres classiques. Les couples $(U_i, \lambda_i)$ sont ici construits les uns après les autres. Les étapes 3 à 7 correspondent à une résolution itérative du système ((4),(5)) pour la recherche d'un nouveau couple $(U_i, \lambda_i)$. Peu d'itérations suffisent généralement à obtenir une bonne approximation de ce couple. En pratique, on choisit $k_{max} = 3$. L'étape 8, qui correspond à une réactualisation des fonctions stochastiques, est facultative. Cependant, elle permet souvent de réduire le nombre $M$ de couples nécessaires pour une précision donnée.

On notera que l'algorithme proposé nécessite d'une part la résolution de problèmes déterministes (étape 4) et d'autre part de problèmes stochastiques sur une base réduite de fonctions déterministes (étapes 6 et 8). Cela conduit à un coût de calcul considérablement réduit comparé aux techniques de résolution classiques.

REMARQUE. — Dans le cas d'un opérateur $\mathbf{A}$ symétrique déterministe, on montre que l'approximation est solution d'un problème aux valeurs propres généralisé classique $\mathbf{A}\mathbf{W} = E(\mathbf{b}\mathbf{b}^T)\mathbf{W}\mathbf{D}$ où $\mathbf{D}$ est la
matrice diagonale des valeurs propres. La meilleure décomposition spectrale est alors une décomposition spectrale classique du second membre b dans la métrique induite par A. L'algorithme 1 correspond dans ce cas à une méthode des puissances classique [Saad (1992)] pour la résolution du problème aux valeurs propres. On notera que dans ce cas, la phase 8 de réactualisation est sans effet.

**Algorithmme 1 Power-type Generalized Spectral Decomposition with Updating (PU-GSD)**

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

2: Initialisation $\lambda_i \in S_P$

3: **for** k = 1 ... $k_{max}$ **do**

4: Calcul de $\mathbf{U}_i$ (problème déterministe) : $\mathbf{U}_i := E(\mathbf{A}\boldsymbol{\lambda}^2)^{-1}E(\mathbf{b}\boldsymbol{\lambda})$

5: Normalisation de $\mathbf{U}_i$ et vérification de la convergence

6: Calcul de $\lambda_i$ (équation stochastique) : $E(\tilde{\lambda}(\mathbf{U}_i^T \mathbf{A} \mathbf{U}_i)\lambda_i) = E(\tilde{\lambda}(\mathbf{U}_i^T \mathbf{b}))$, $\forall \tilde{\lambda} \in S_P$

7: **end for**

8: Poser $\mathbf{W} = (\mathbf{U}_1 \dots \mathbf{U}_i)$ et réactualisation de $\mathbf{\Lambda}$ en résolvant (5) (système réduit d'équations stochastiques)

9: Réactualisation de $\mathbf{b}$ et vérification de la convergence

10: **end for**

**4 Problèmes non-linéaires : réutilisation des bases spectrales**

Dans le cas de problèmes non-linéaires, la formulation discrétisée du problème peut s'écrire :

$$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 (6)$$

où b est un vecteur aléatoire et où $\mathcal{A}$ est un opérateur non-linéaire. Ce système est classiquement résolu par un algorithme itératif. En notant $\mathbf{u}^{(k)}$ l'approximation à l'itération k, l'itéré suivant $\mathbf{u}^{(k+1)} \in \mathbb{R}^n \otimes S_P$ est solution du système :

$$E(\mathbf{v}^T \hat{\mathbf{A}}^{(k)} \mathbf{u}^{(k+1)}) = E(\mathbf{v}^T (\mathbf{b} - \mathcal{A}(\mathbf{u}^{(k)})) + \hat{\mathbf{A}}^{(k)} \mathbf{u}^{(k)}) \quad \forall \mathbf{v} \in \mathbb{R}^n \otimes S_P \quad (7)$$

où $\hat{\mathbf{A}}^{(k)}$ est une matrice dépendant du solveur non-linéaire utilisé. Dans le cas d'un algorithme de Newton-Raphson classique, cette matrice est classiquement prise comme la dérivée de Gâteaux de $\mathcal{A}$, définie $\forall \mathbf{w} \in \mathbb{R}^n \otimes S_P$ par : $\hat{\mathbf{A}}^{(k)}\mathbf{w} = \lim_{\epsilon \to 0} \frac{1}{\epsilon}(\mathcal{A}(\mathbf{u}^{(k)} + \epsilon\mathbf{w}) - \mathcal{A}(\mathbf{u}^{(k)}))$. La résolution du problème non-linéaire se ramène donc à la résolution d'un ensemble de problèmes stochastiques linéaires du type (1). L'algorithme 1 peut ainsi être utilisé pour la résolution de chacun de ces systèmes. Un grand intérêt de l'approche présentée est que les bases spectrales généralisées construites au cours des itérations peuvent être stockées et réutilisées pour les résolutions succéssives. Ces bases s'avèrent très pertinentes vis à vis de l'opérateur et du second membre du problème. La réutilisation devient donc très efficace lorsque l'opérateur et le second membre des problèmes linéaires successifs varient peu d'une itération sur l'autre.

**5 Exemple : problème de thermique non-linéaire**

On considère un problème de thermique non-linéaire défini sur le domaine $\Omega$ (figure 1). Des conditions homogènes de Dirichlet sont imposées sur la frontière $\partial\Omega$. On utilise au niveau spatial une approximation élément finis $\mathcal{P}_1$ (quadrangle à 4 noeuds) avec $n = 1121$ degrés de liberté. La température $u(x, \theta)$ semi-discrétisée au niveau déterministe est recherchée sous la forme $u(x, \theta) = \varphi(x)u(\theta)$, $\mathbf{u} \in \mathbb{R}^n \otimes S_P$, où $\varphi$ désigne le vecteur ligne des fonctions d'interpolation éléments finis. L'opérateur et le second membre
du problème (6) s'écrivent

$$E(\mathbf{v}^T \mathbf{A}(\mathbf{u})) = E(a(\varphi \mathbf{u}, \varphi \mathbf{v}; \theta)) \quad \text{avec} \quad a(u, v; \theta) = \int_{\Omega} \nabla v \cdot (\kappa_0(\theta) + \kappa_1(\theta) u^2) \nabla u \, dx$$

$$E(\mathbf{v}^T \mathbf{b}) = E(b(\varphi \mathbf{v}; \theta)) \quad \text{avec} \quad b(v; \theta) = \int_{\Omega} f(\theta) v \, dx$$

où $f$ est un terme source de chaleur homogène et $\kappa_0$ et $\kappa_1$ des paramètres matériau. On modélise $f$, $\kappa_0$ et $\kappa_1$ par 3 variables aléatoires indépendantes : $f(\theta) = \mathcal{N}(1, 0.2)^1$, $\kappa_0 = U(0.7, 1.3)$ et $\kappa_1 = U(0.7, 1.3)^2$. On choisit au niveau stochastique une approximation par chaos polynômial généralisé de degré $p = 3$ en dimension $m = 3$ [Soize et al. (2004)]. L'espace d'approximation $S_P$ est donc de dimension $P = \frac{(m+p)!}{m!p!}$.

Le problème est résolu par un algorithme de Newton-Raphson avec une tolérance de $10^{-4}$. Chaque système linéaire (7) est résolu avec une précision relative de $10^{-5}$. L'algorithme (PU-GSD) proposé est comparé avec un algorithme classique Conjugate Gradient Square (CGS) [Saad (1992)]. Dans le cas où les bases spectrales sont réutilisées, on applique systématiquement une procédure de tri basée sur une décomposition spectrale classique.

<figure><img src="image_2.png" /><figcaption>FIG. 1 – Présentation du problème</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 2 – Modes {$U_i$}$_{i=1}^6$ obtenues par PU-GSD avec réutilisation</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 3 – Nombres de vecteurs déterministes $U_i$ construits par PU-GSD au cours des itérations de Newton-Raphson : influence de la réutilisation</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4 – Temps de calcul pour la première itération de Newton-Raphson en fonction de la dimension $P$ de l'espace d'approximation stochastique</figcaption></figure>

L'algorithme (PU-GSD), avec réutilisation et tri des modes déterministes, conduit à convergence à une décomposition spectrale en $M = 10$ modes. La figure 2 montre les 6 premiers modes. On notera

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

<footnote n="2">²$U(\alpha, \beta)$ désigne une variable uniforme sur $[\alpha, \beta]$.</footnote>
que ces 6 modes sont les plus significatifs puisque que la solution spectrale tronquée à $M = 6$ conduit à une erreur de $10^{-4}$. La figure 3 montre le nombre de couples $(U_i, \lambda_i)$ construits à chaque itération de Newton. On remarque que pour la résolution du premier système linéaire, un seul couple suffit à représenter la solution exacte, ce qui était prévisible. On remarque également qu'avec la réutilisation des bases spectrales, de moins en moins de nouveaux couples sont nécessaires au fil des itérations de Newton. Enfin, la figure 4 montre l'influence du degré du chaos polynômial $p$, et donc de la dimension stochastique $P$, sur le temps de calcul de l'ensemble des systèmes linéaires. On remarque que l'algorithme proposé est très peu sensible à la dimension stochastique tandis que le temps de calcul d'une approche standard (CGS) croît très fortement avec la dimension $P$.

# 6 Conclusion

Nous avons proposé une méthode de résolution de problèmes stochastiques qui consiste à générer automatiquement une base spectrale généralisée de la solution. Les avantages de cette approche sont multiples : réduction drastique du coût de stockage mémoire de la solution, réduction du temps de calcul par découplage de problèmes déterministes et de problèmes stochastiques sur bases réduites, possibilité de réutiliser efficacement les bases spectrales généralisées au cours d'une résolution itérative d'un problème non-linéaire ou plus généralement dans le cadre de multi-résolutions.

## Remerciements

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

# Références

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

M. Deb, I. Babuška, and J. T. Oden. Solution of stochastic partial differential equations using galerkin finite element techniques. *Computer Methods in Applied Mechanics and Engineering*, 190 :6359–6372, 2001.

H. G. Matthies and A. Keese. Galerkin methods for linear and nonlinear elliptic stochastic partial differential equations. *Comp. Meth. App. Mech. Eng.*, 194(12-16) :1295–1331, 2005.

R. G. Ghanem and R. M. Kruger. Numerical solution of spectral stochastic finite element systems. *Computer Methods in Applied Mechanics and Engineering*, 129 :289–303, 1996.

M. F. Pellissetti and R. G. Ghanem. Iterative solution of systems of linear equations arising in the context of stochastic finite elements. *Advances in Engineering Software*, 31 :607–616, 2000.

A. Keese and H. G. Matthies. Hierarchical parallelisation for the solution of stochastic finite element equations. *Comp. Meth. App. Mech. Eng.*, 83 :1033–1047, 2005.

A. Nouy. A generalized spectral decomposition technique to solve stochastic partial differential equations. *Submitted to Computer Methods in Applied Methods in Engineering*, 2007.

C. Soize and R. Ghanem. Physical systems with random uncertainties : chaos representations with arbitrary probability measure. *SIAM J. Sci. Comput.*, 26(2) :395–410, 2004.

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

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