# Décomposition en série de Karhunen-Loève d'un processus aléatoire par le biais d'un code éléments finis

Sébastien Recek, Maurice Lemaire, Alain Millard

► To cite this version:

Sébastien Recek, Maurice Lemaire, Alain Millard. Décomposition en série de Karhunen-Loève d'un processus aléatoire par le biais d'un code éléments finis. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01813022

HAL Id: hal-01813022

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

Submitted on 12 Jun 2018

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
# Décomposition en Série de Karhunen-Loève d'un Processus Aléatoire par le Biais d'un Code Eléments Finis

Sébastien Recek* — Maurice Lemaire** — Alain Millard *

* Commissariat à l'Energie Atomique
DM2S / SEMT / LM2S
CEA / SACLAY
F-91191 Gif-sur-Yvette

** LaMI - UBP & IFMA
EA 3867 / FR CNRS 2856
F-63175 Aubière

RÉSUMÉ. Ce papier présente l'implémentation d'une décomposition numérique de Karhunen-Loève d'une fonction de covariance dans un code de calcul par éléments finis. Cette procédure permet la discrétisation future d'un champ aléatoire sous la forme d'une série tronquée comportant des variables aléatoires et des fonctions déterministes.

ABSTRACT. This paper presents the implementation of a numerical Karhunen-Loève covariance expansion method in a finite elements code. The procedure allows the future discretization of a stochastic field in the form of a truncated series expansion including deterministic functions and random variables.

MOTS-CLÉS : Eléments Finis Stochastiques, Champ Aléatoire, Processus Aléatoire, Décomposition Spectrale, Décomposition en série de Karhunen-Loève, Covariance.

KEYWORDS: Stochastic Finite Elements, Stochastic fields, Stochastic processes, Spectral expansion, Karhunen-Loève Expansion, Covariance.
# 1. Introduction

L'une des difficultés principales associée aux calculs stochastiques est la manipulation de champs aléatoires, utilisés par exemple pour représenter la variabilité spatiale de données (densité, module d'Young,...). Il est nécessaire de les discrétiser, les approximer sous une forme comprenant un nombre fini de variables aléatoires. Pour un champ aléatoire cela peut par exemple consister en l'association d'une variable aléatoire à chaque nœud d'un modèle éléments finis et en l'utilisation des fonctions de forme pour interpoler le champ aléatoire sur tout le domaine. Cette méthode peut conduire à un grand nombre de variables aléatoires, avec lequel il sera coûteux de travailler (pour réaliser des simulations de Monte Carlo ou des calculs stochastiques).

Une démarche moins coûteuse consiste en la décomposition du champ sous la forme d'une série analogue à une série de Fourier [GHA 91]. Considérons un champ aléatoire $\alpha(\{x_i\}, \omega)$ dont sont connues la moyenne $\bar{\alpha}(\{x\})$ et la fonction de covariance $C_\alpha(\{x\}, \{y\})$. La décomposition en série de Karhunen-Loève du processus aléatoire est de la forme :

$$ \alpha(\{x\}, \omega) = \bar{\alpha}(\{x\}) + \sum_{n=1}^{\infty} \sqrt{\lambda_n} \zeta_n(\omega) \phi_n(\{x\}) \quad [1] $$

Où $\lambda_n$ et $\phi_n(\{x\})$ sont les valeurs propres et vecteurs propres de la fonction de covariance, c'est à dire les solution de l'équation de Fredholm suivante (D étant le domaine d'intégration et $d\{x_i\}$ devant être compris comme le produit des $dx_{ij}$) :

$$ \int_D C_\alpha(\{x_1\}, \{x_2\}) \phi_n(\{x_1\}) d\{x_1\} = \lambda_n \phi_n(\{x_2\}) \quad [2] $$

Les variables aléatoires $\zeta_n(\omega)$ Normales si le champ aléatoire $\alpha(\{x\}, \omega)$ est Gaussien. La convergence de la série [1] dépend de la fonction de covariance [HUA 01]. Une telle décomposition peut constituer la première étape d'un calcul par une méthode d'éléments finis stochastiques dans le cadre de processus aléatoire Gaussien ([GHA 99], [MAT 99]).

Des solutions analytiques à l'équation [2] sont connues pour quelques fonctions de covariances classiques et des domaines unidimensionnels. Cependant dans [GHA 91] une méthode de résolution numérique est présentée permettant la résolution de [2] pour tout domaine et fonction de covariance admissible. Cette méthode présente l'avantage d'être particulièrement adaptée à une implémentation dans un code de calcul par éléments finis, même si elle peut être utilisée hors de ce cadre [CAS 02]. Dans ce papier nous présentons une implémentation de cette méthode dans le code de calcul par éléments finis CAST3M [VER 84] et la testons sur un cas possédant une solution analytique.
## 2. Théorie

Une procédure de type Galerkin permet d'obtenir une solution à l'équation [2] [GHA 91]. Soit $\{h_i(\{x\})\}$ une base de fonctions dans l'espace de Hilbert des fonctions continues. Chaque fonction propre $\phi_k(\{x\})$ se décompose en :

$$ \phi_k(\{x\}) = \sum_{i=1}^{\infty} d_i^{(k)} h_i(\{x\}) \quad [3] $$

La différence $\varepsilon_k$ entre les deux membres de l'équation de Fredholm [2] suite à la troncature de cette série au $N^{\text{ième}}$ terme pour la $k^{\text{ième}}$ fonction propre est :

$$ \varepsilon_k(\{x_1\}) = \sum_{i=1}^{N} d_i^{(k)} \left[ \int_D C_\alpha(\{x_1\}, \{x_2\}) h_i(\{x_2\}) d\{x_2\} - \lambda_k h_i(\{x_1\}) \right] \quad [4] $$

Pour un $k$ donné requérir que l'erreur soit orthogonale aux fonctions de la base conduit à :

$$ \sum_{i=1}^{N} d_i^{(k)} \left[ \int_D \left( \int_D C_\alpha(\{x_1\}, \{x_2\}) h_i(\{x_2\}) d\{x_2\} \right) h_i(\{x_1\}) d\{x_1\} - \lambda_k \int_D h_i(\{x_1\}) h_i(\{x_1\}) d\{x_1\} \right] = 0 \quad [5] $$

Cette équation permet de calculer la $k^{\text{ième}}$ fonction et valeur propre. Cette expression peut être écrite sous une forme matricielle en définissant les matrices [C] et [B] suivantes :

$$ C_{ij} = \int_D \left( \int_D C_\alpha(\{x_1\}, \{x_2\}) h_i(\{x_2\}) d\{x_2\} \right) h_j(\{x_1\}) d\{x_1\} \quad [6] $$

$$ B_{ij} = \int_D h_i(\{x\}) h_j(\{x\}) d\{x\} \quad [7] $$

Ce qui conduit à l'équation suivante pour le $k^{\text{ième}}$ mode :

$$ [B]\{d^{(k)}\} - \frac{1}{\lambda_k} [C]\{d^{(k)}\} = 0 \quad [8] $$

Cette équation est similaire à celle définissant un problème dynamique dans le cadre d'une méthode de calcul par éléments finis. En effet, pour un système possédant une matrice de rigidité [K] et de masse [M], les $k^{\text{ième}}$ pulsation et vecteur propre sont obtenu en résolvant l'équation :

$$ [K]\{q^{(k)}\} - \omega_k^2 [M]\{q^{(k)}\} = 0 \quad [9] $$

Par conséquent il suffit d'ajouter au code le calcul des matrices [C] et [B] pour disposer d'un outil de résolution de l'équation [2]. Notre procédure de calcul conduit aux correspondances suivantes dans le modèle :
Données :

- la matrice [K] est remplacée par la matrice [B] :

- la matrice [M] est remplacée par la matrice [C] ;

- Aucune condition aux limites n'est imposée au système

Résultats:

- la k-ième pulsation $\omega_k$ retournée nous donne la valeur propre $\lambda_k$ selon la relation :

$$ \lambda_k = \frac{1}{\omega_k^2} \qquad [10] $$

- le vecteur propre {$q^{(k)}$} retourné suite au calcul nous donne les valeurs nodales de la fonction propre $\phi_k(\{x\})$ qui peut ensuite être interpolée en utilisant les n fonctions de forme $N_i(\{x\})$ (n étant le nombre de nœuds du modèle) :

$$ \phi_k(\{x\}) = \sum_{i=1}^{n} N_i(\{x\}) q_i^{(k)} \qquad [11] $$

### 3.Calcul des matrices [C] et [B] en utilisant un modèle éléments finis

L'utilisation des fonctions de forme $N_i(\{x\})$ comme base de fonction dans l'expression [6] conduit à :

$$ C_{ij} = \int_D \left[ \int_D C(\{x_1\}, \{x_2\}) N_i(\{x_2\}) d\{x_2\} \right] N_j(\{x_1\}) d\{x_1\} \qquad [12] $$

Considérons d'abord le premier niveau d'intégration de l'équation [6] :

$$ f_i(\{x_1\}) = \int_D C(\{x_1\}, \{x_2\}) N_i(\{x_2\}) d\{x_2\} \qquad [13] $$

Les fonctions de forme possédant un support compact, le domaine d'intégration peut être limité aux éléments entourant directement le nœud i. En notant Di ce domaine et De_k le domaine de l'élément k nous permet d'écrire :

$$ f_i(\{x_1\}) = \sum_{k_1 \in D_i} \left( \int_{D_{k_1}} C(\{x_1\}, \{x_2\}) N_i(\{x_2\}) d\{x_2\} \right) \qquad [14] $$

Puis, en utilisant une intégration par point de Gauss :

$$ f_i(\{x_1\}) = \sum_{k_1 \in D_i} \sum_{ig=1}^{n_{\text{gauss}}} C(\{x_1\}, \{x_{ig}^{k_1}\}) N_i(\{x_{ig}^{k_1}\}) w_{ig} \det([J_{k_1}]) \qquad [15] $$

De même, pour la variable {$x_1$} :

$$ C_{ij} = \int_{D_i} f_i(\{x_1\}) N_j(\{x_1\}) d\{x_1\} \qquad [16] $$
$$C_{ij} = \sum_{k_2 \in D_i} \sum_{jg=1}^{n_{\text{Gauss}}} f(\{x_{jg}^{k_2}\}) N_j(\{x_{jg}^{k_2}\}) w_{jg} \det([J_{k_2}]) \quad [17]$$

Finalement, en substituant f par son expression [15], nous obtenons :

$$C_{ij} = \sum_{k_2 \in D_i} \sum_{jg=1}^{n_{\text{Gauss}}} \left[ \sum_{k_1 \in D_i} \sum_{ig=1}^{n_{\text{Gauss}}} C(\{x_{jg}^{k_2}\}, \{x_{ig}^{k_1}\}) N_i(\{x_{ig}^{k_1}\}) w_{ig} \det([J_{k_1}]) \right] N_j(\{x_{jg}^{k_2}\}) w_{jg} \det([J_{k_2}]) \quad [18]$$

Où Di est le jeu d'éléments entourant le nœud i, $N_i$ est la fonction de forme associée au nœud i, $n_{\text{Gauss}}$ est le nombre de points de Gauss dans chaque élément, $\{x_{ig}^{k_1}\}$ est le vecteur des coordonnées du $ig^{th}$ point de Gauss de l'élément $k_1$, $w_{ig}$ est le poids associé au $ig^{th}$ point de Gauss $[J_k]$ est la matrice Jacobienne de l'élément $k_1$.

Pour la matrice [B] le calcul est plus simple : l'expression [7] correspond à la matrice de masse du système pour une densité unitaire.

## 6. Cas-test

Nous avons utilisé comme cas-test une poutre maillée avec des éléments tétraédriques et une covariance exponentielle, pour lesquels une solution analytique est connue. Les résultats numériques montrent une bonne convergence vers les résultats analytiques avec l'augmentation de la densité du maillage (Figure 1). Pour la densité la plus élevée testée, les valeurs propres et fonctions propres (Figure 2) numériques pour les 4 premiers modes sont similaires à 1% près à la théorie.

<figure><img src="image_2.png" /><figcaption>Figure 1. Cas-test : erreur relative entre les valeurs propres calculées et les valeurs propres théoriques en fonction du nombre d'éléments</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Cas-test : modes propres calculés (à gauche) et théoriques (à droite)</figcaption></figure>

**8. Conclusions**

La méthode présentée précédemment permet la décomposition en série de Karhunen-Loève d'une fonction de covariance par le biais d'un code de calcul par éléments finis. Cette décomposition permet la discrétisation de champs aléatoires sous la forme d'un jeux de fonctions déterministes (représentées par leurs valeurs nodales) associées à des variables aléatoires. Cette méthode est utilisable pour toute fonction de covariance admissible et tout domaine.

**References**

[GHA 91] Ghanem R, Spanos P. “Stochastic Finite Elements: A Spectral Approach”. Springer-Verlag 1991.

[HUA 01] Huang SP, Quek ST, Phoon KK. “Convergence Study of the Truncated Karhunen-Loève Expansion for simulation of Stochastic Processes.” *Int. J. Numer. Meth. Engng* 2001; 52:1029-1043.

[GHA 99] Ghanem R. “Ingredients for a General Purpose Stochastic Finite Elements Implementation”. *Comp. Methods Appl. Engrg.* 1999; 168:19-34.

[MAT 99] Matthias HG, Bucher C. “Finite Elements for stochastic Media Problems”. *Comp. Methods Appl. Engrg.* 1999; 168:19-34.

[CAS 02] Castrillon-Candas JE, Amaratunga K. “Fast Estimation of continuous Karhunen-Loève Eigenfunctions Using Wavelets.” *IEEE Transaction on signal Processing* 2002; 50(1):78-86.

[VER 84] Verpeaux P, Millard A, Charras T, Combescure A. “A Modern Approach of Large Computer Codes for Structural Analysis”. *Proc. SMIRT 10*, Los Angeles, 1984.