ANNALES MATHÉMATIQUES BLAISE PASCAL

PH. BRESSOLETTE

M. FOGLI

Quelques méthodes de résolution de problèmes de dynamique stochastique non linéaire

*Annales mathématiques Blaise Pascal*, tome 1, n° 1 (1994), p. 43-60

<http://www.numdam.org/item?id=AMBP_1994__1_1_43_0>

© Annales mathématiques Blaise Pascal, 1994, tous droits réservés.

L'accès aux archives de la revue « Annales mathématiques Blaise Pascal » (http://math.univ-bpclermont.fr/ambp/) implique l'accord avec les conditions générales d'utilisation (http://www.numdam.org/legal.php). Toute utilisation commerciale ou impression systématique est constitutive d'une infraction pénale. Toute copie ou impression de ce fichier doit contenir la présente mention de copyright.

NUMDAM

Article numérisé dans le cadre du programme

Numérisation de documents anciens mathématiques

http://www.numdam.org/
# QUELQUES MÉTHODES DE RÉSOLUTION DE PROBLÈMES DE DYNAMIQUE STOCHASTIQUE NON LINÉAIRE

Ph. BRESSOLETTE - M. FOGLI

Cet article est issu du séminaire Fiabilité, tenu le 29 janvier 1992 à Clermont-Ferrand, dans le cadre d'une journée d'étude de la Formation Doctorale MATÉRIAUX, STRUCTURES, FIABILITÉ.

## 1. INTRODUCTION

Nous nous proposons dans cet article d'exposer quelques méthodes de résolution pour les problèmes de mécanique des vibrations aléatoires non linéaires. Ces problèmes se rendent classiquement en fiabiblité structurale (étude des effets aléatoires du vent, de la houle ou des séismes sur les structures, fatigue sous chargements aléatoires, modélisation des mouvements des milieux granulaires, vibrations d'éléments de transport de fluides turbulents...). Ils se ramènent généralement à la résolution de systèmes dynamiques vectoriels de dimension finie, non linéaires, du second ordre et à excitation aléatoire externe et paramétrique, de la forme :

$$ egin{cases} \ddot{Q}(t) + s(Q(t), \dot{Q}(t), E(t), \Lambda(t)) = 0, & t > 0 \\ Q(0) = Q_0 \text{ p.s. ; } \dot{Q}(0) = \dot{Q}_0 \text{ p.s.} \end{cases} \tag{1} $$

où s est une application de $\mathbb{R}^n \times \mathbb{R}^n \times \mathbb{R}^n \times \mathbb{R}^p$ dans $\mathbb{R}^n$, $E$ est le processus excitation externe à valeurs dans $\mathbb{R}^n$, $\Lambda$ le processus excitation paramétrique à valeurs dans $\mathbb{R}^p$, indépendante de $E$, $Q$ le processus réponse à valeurs dans $\mathbb{R}^n$ et $X_0 = (Q_0, \dot{Q}_0)^t$ une variable aléatoire (v.a.) à valeurs dans $\mathbb{R}^n \times \mathbb{R}^n = \mathbb{R}^{2n}$, indépendante de $\{E(t), t \ge 0\}$ et $\{\Lambda(t), t \ge 0\}$, et de loi $P_{X_0}$ donnée sur $(\mathbb{R}^{2n}, \mathcal{B}_{2n})$. Le processus $Q$ représente par exemple le vecteur des déplacements nodaux si (1) décrit les vibrations aléatoires d'un oscillateur à n degrés de liberté (ddl).

Afin de simplifier l'exposé, nous ne considérerons ici que des systèmes scalaires (i.e. $n = 1$) à excitation externe. Nous supposerons en outre que l'excitation est gaussienne large bande, donc modélisable par un bruit blanc gaussien, et que la loi de la donnée initiale admet une densité par rapport à la mesure de Lebesgue. De tels systèmes peuvent alors s'écrire sous la forme :

$$ egin{cases} \ddot{Q}(t) + f(Q(t), \dot{Q}(t)) = D\dot{W}(t), & t > 0 \\ Q(0) = Q_0 \text{ p.s. ; } \dot{Q}(0) = \dot{Q}_0 \text{ p.s.} \end{cases} \tag{2} $$
où $W$ est le processus de Wiener normalisé à valeurs dans $\mathbb{R}$, $\dot{W}$ sa dérivée distribution par rapport au temps, c'est-à-dire le bruit blanc gaussien normalisé scalaire (de forme de covariance $E[\dot{W}(t+\tau)\dot{W}(t)] = \delta_0(\tau)$ et de densité spectrale $S_{\dot{W}}(\omega) = 1/2\pi; \omega \in \mathbb{R}$), $f$ est une application de $\mathbb{R}^2$ dans $\mathbb{R}$, $D$ un réel strictement positif et $X_0 = (Q_0, \dot{Q}_0)^t$ une v.a. à valeurs dans $\mathbb{R}^2$, indépendante de $\{W(t), t \ge 0\}$ et de loi $P_{X_0}$ admettant une densité $p_{X_0}$ par rapport à la mesure de Lebesgue sur $\mathbb{R}^2: P_{X_0}(dx_0) = p_{X_0}(x_0)dx_0$; $x_0 = (q_0, \dot{q}_0)$, $dx_0 = dq_0d\dot{q}_0$.

Notons qu'en passant dans l'espace des phases, c'est-à-dire en introduisant le processus $X$ tel que $X = (Q, \dot{Q})^t$, l'équation (2) se ramène classiquement à une équation différentielle stochastique d'Itô (EDSI) sur $\mathbb{R}^2$ de la forme :

$$ \begin{cases} dX(t) = b(X(t)) dt + \sigma(X(t)) dW(t), & t > 0 \\ X(0) = X_0 \quad \text{p.s.} \end{cases} \qquad (3) $$

avec, pour tout $x = (q, \dot{q})$ dans $\mathbb{R}^2$ :

$$ b(x) = \begin{pmatrix} \dot{q} \\ -f(q, \dot{q}) \end{pmatrix} ; \sigma(x) = \sigma = \begin{pmatrix} 0 \\ D \end{pmatrix} \qquad (4) $$

Par suite, résoudre (2) équivaut à résoudre (3). Signalons au passage que le fait de considérer une excitation blanche n'est pas une perte de généralité. En effet, si l'excitation était gaussienne colorée et physiquement réalisable, on pourrait toujours se ramener à un problème de type (3), mais de dimension supérieure, en utilisant une technique de blanchiment bien connue : la méthode de markovianisation [6,16].

Dans tout ce qui suit, nous nous intéressons essentiellement à la solution stationnaire de (2)-(3), dont nous admettrons l'existence et l'unicité et que nous supposerons du second ordre. Cela implique que (2.1) ou (3.1) admet une unique probabilité invariante du second ordre et que le problème d'évolution (2) ou (3) admet une unique solution qui n'explose pas. Nous supposerons en outre que $f$ vérifie $f(-q, -\dot{q}) = -f(q, \dot{q})$, propriété qui assurera à la solution d'être centrée.

En fait, dès que la dimension de l'espace des phases (i.e. la dimension de $X$) dépasse 1, c'est-à-dire dans tous les cas pratiques, on est incapable de résoudre exactement et par voie purement analytique les problèmes associés à la construction des solutions des EDSI non linéaires (d'où, notamment pour les réponses stationnaires, l'impossibilité d'obtenir les expressions analytiques exactes des densités spectrales). Or ces quantités sont d'intérêt pour le mécanicien. Elles seules, en effet, permettent de montrer comment se répartit l'énergie de la réponse par bande de fréquence. En outre, l'étude de leur comportement asymptotique donne de précieuses indications sur la régularité du processus réponse. Il importe donc, à défaut de solutions explicites exactes, de disposer au moins d'approximations.

Les méthodes que nous proposons répondent à cette exigence et sont, du moins en principe pour certaines d'entre elles, généralisables à la dimension supérieure. Bien entendu, nous ne prétendons pas à l'exhaustivité, non plus qu'à l'objectivité dans le choix des méthodes. Nous présenterons successivement les méthodes basées sur la simulation
numérique, la méthode des processus de diffusion, la méthode des perturbations et les méthodes de linéarisation stochastique.

Concernant l'exposé, nous nous bornerons à un bref rappel du principe de chaque méthode, laissant de côté les justifications théoriques et les développements constructifs qui leurs sont attachés. Nous proposons toutefois en fin d'article une bibliographie sélectionnée dans laquelle le lecteur intéressé pourra trouver tous ces détails.

## 2. LES MÉTHODES BASÉES SUR LA SIMULATION NUMÉRIQUE ([1] à [9])

Ces méthodes consistent à construire artificiellement, par simulation numérique de type Monte-Carlo, des trajectoires de la solution stationnaire de (3), puis à traiter statiquement ces trajectoires à l'aide de méthodes de traitement du signal afin d'en déduire des estimations des densités spectrales. La simulation numérique des trajectoires de la diffusion stationnaire solution de l'EDSI (3) nécessite la mise en œuvre d'algorithmes qui s'obtiennent en discrétisant (3) à l'aide de schémas numériques appropriés.

Divers types de schémas peuvent être utilisés. Par exemple, en utilisant un schéma implicite aux différences centrées de pas $\Delta$ constant, la discrétisation de (3) donne, en posant $X(t_i) = X_i, t_i = i\Delta, i = 0, N$ :

$$ \begin{cases} X_0 \\ X_{i+1} = X_i + b \left( \frac{X_{i+1} + X_i}{2} \right) \Delta + \sqrt{\Delta} \sigma G_i; & i = 0, N \end{cases} \quad (5) $$

où $(G_i; i = 0, N)$ est une suite de copies indépendantes d'une v.a. $G$ gaussienne réduite à valeurs dans $\mathbb{R}$ et $X_0$ une v.a. à valeurs dans $\mathbb{R}^2$, indépendante de $G$, et ayant pour loi la probabilité invariante de (3.1) (cf § 3).

Par suite, l'algorithme s'écrit :

$$ \begin{cases} x_0 \\ x_{i+1} = F_{x_i, g_i}(x_{i+1}); & i = 0, N \end{cases} \quad (6) $$

où $x \to F_{x_i,g_i}(x)$ est l'application de $\mathbb{R}^2$ dans $\mathbb{R}^2$ définie par :

$$ F_{x_i, g_i}(x) = x_i + b \left( \frac{x+x_i}{2} \right) \Delta + \sqrt{\Delta} \sigma g_i \quad (7) $$

($g_i; i = 1, N$) étant une suite de réalisations indépendantes de $G$ et $x_0$ une réalisation de $X_0$.

Avec cet algorithme, on est donc ramené, à chaque instant de la discrétisation, à la résolution d'un problème de point fixe sur $\mathbb{R}^2$ du type : $x = F_{x_i,g_i}(x)$.

On pourrait également discrétiser (3) à l'aide d'un schéma explicite d'Euler. On obtiendrait alors, en conservant les mêmes notations que précédemment :

$$ \begin{cases} X_0 \\ X_{i+1} = X_i + b(X_i)\Delta + \sqrt{\Delta}\sigma G_i; & i = 0, N \end{cases} \quad (8) $$
et donc pour l'algorithmme associé :

$$ \begin{cases} x_0 \\ x_{i+1} = F_{g_i}(x_i) ; i = 0, N \end{cases} \qquad (9) $$

avec $x \to F_{g_i}(x)$ l'application de $\mathbb{R}^2$ dans $\mathbb{R}^2$ définie pour tout $x$ par :

$$ F_{g_i}(x) = x + b(x)\Delta + \sqrt{\Delta} \sigma g_i \qquad (10) $$

Toutefois, contrairement aux algorithmes implicites, les algorithmes explicites de ce type présentent l'inconvénient de nécessiter de très petits pas de discrétisation (problèmes de stabilité numérique) et par conséquent d'élever le coût des simulations. Bien que conduisant à chaque étape de la simulation à la résolution d'une équation algébrique, il est donc préférable de choisir des schémas numériques implicites (qui ne permettent cependant pas, pour des raisons de précision, un choix quelconque du pas de discrétisation).

L'estimation des densités spectrales se fait quant à elle classiquement à partir des trajectoires simulées de la réponse en utilisant la méthode du périodogramme avec fenêtre temporelle [3,6].

### 3. LA MÉTHODE DES PROCESSUS DE DIFFUSION ([10] à [22])

Cette méthode opère sur la formulation (3) du problème, qui est une équation de diffusion de vecteur de dérive $b(x) \in \mathbb{R}^2$ et de matrice de diffusion $a(x) = \sigma(x)\sigma(x)^t \in Mat_R(2, 2)$. Comme les fonctions $a$ et $b$ ne dépendent pas explicitement du temps, le processus de Markov $X$ solution de (3) est homogène. Il est complètement caractérisé par la donnée de la loi initiale $P_{X_0}$ et par sa probabilité de transition $P(X(t) \in B/X(0) = x_0)$, $B$ borélien de $\mathbb{R}^2$, $x_0 \in \mathbb{R}^2$, dont nous supposerons qu'elle admet une densité $\rho(x, t/x_0)$. Dans ces conditions, $\rho$ vérifie l'équation de Fokker-Planck instationnaire (EFPI) :

$$ \begin{cases} \frac{\partial \rho}{\partial t} - L^T \rho = 0, & t > 0 \\ \rho(x, t/x_0) dx \to \delta_{x_0}(dx) & \text{si } t \downarrow 0 \end{cases} \qquad (11) $$

où $\delta_{x_0}$ est la mesure de Dirac au point $x_0$ et $L^T$ l'opérateur linéaire défini par :

$$ L^T \rho = - \sum_{i=1}^{2} \frac{\partial}{\partial x_i} (b_i(x) \rho) + \frac{1}{2} \sum_{i=1}^{2} \sum_{j=1}^{2} \frac{\partial^2}{\partial x_i \partial x_j} (a_{ij}(x) \rho) \qquad (12) $$

De plus, d'après les hypothèses, la solution de (3) tend asymptotiquement en probabilité pour $t \to +\infty$ vers un processus stationnaire $X_S$ du second ordre à valeurs dans $\mathbb{R}^2$. Par suite, la loi $P_{X_S(t)}$ de la v.a. $X_S(t)$ est indépendante de $t$ et si elle admet une densité $p_S$, ce que nous supposerons, alors celle-ci est solution de l'équation de Fokker-Planck stationnaire (EFPS) :

$$ \begin{cases} L^T p_S = 0 \\ \int_{\mathbb{R}^2} p_S(x) dx = 1 \end{cases} \qquad (13) $$
et vérifie en outre la relation d'invariance :

$$ \forall t \ge 0, \quad p_S(x) = \int_{\mathbb{R}^2} \rho(x, t/x_0) p_S(x_0) dx_0 \quad (14) $$

$P_{X_S(t)}$ est la probabilité invariante de (3.1). Si l'on prend pour condition initiale $X_0$ une v.a. ayant pour loi $P_{X_S(t)}$, c'est-à-dire si l'on prend $p_{X_0}(x_0) = p_S(x_0)$, alors le processus $X_S$ vérifie (3). Donc si l'on s'intéresse uniquement à la solution stationnaire de (3), il faut choisir $X_0$ distribuée suivant la loi invariante de (3.1).

La donnée des densités $\rho$ et $p_S$ caractérise complètement le système des lois finidimensionnelles du processus $X_S$. Rappelons que le système des lois finidimensionnelles de $X_S$ est la famille des lois $P_{X_S^N}$ des v.a. $X_S^N = (X_S(t_1), ..., X_S(t_N))$ à valeurs dans $\mathbb{R}^{2N}$ engendrée lorsque $(t_1, ..., t_N)$ parcourt toutes les parties finies non vides et non ordonnées de $\mathbb{R}_+$.

Comme $X_S$ est markovien, il est possible de construire ce système à partir de la seule connaissance de $\rho$ et $p_S$. Ainsi, pour $t_1 > t_2 > ... > t_N$, la densité de la loi $P_{X_S^N}$ de la v.a. $X_S^N$ s'écrit :

$$ p(x_1, t_1; \dots; x_N, t_N) = p_S(x_1) \rho(x_2, t_2 - t_1/x_1) \dots \rho(x_N, t_N - t_{N-1}/x_N) \quad (15) $$

En particulier, pour tout $t \ge 0$, la densité de la loi du couple $(X(0), X(t))$ a pour expression :

$$ p(x_0, 0; x, t) = p_S(x_0) \rho(x, t/x_0) \quad (16) $$

On peut alors en déduire la fonction de corrélation de la solution stationnaire $Q_S$ de (2), qui s'écrit pour $t \ge 0$ :

$$ R_Q(t) = \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} q q_0 p_S(q_0, \dot{q}_0) \rho(q, \dot{q}, t/q_0, \dot{q}_0) dq_0 d\dot{q}_0 dq d\dot{q} \quad (17) $$

et qui vérifie pour tout $t \in \mathbb{R} : R_Q(t) = R_Q(-t)$.

Par suite, si $R_Q \in L^1(\mathbb{R}, \mathbb{R})$, la densité spectrale de $Q_S$ s'obtient par :

$$ S_Q(\omega) = \frac{1}{2\pi} \int_{\mathbb{R}} e^{-i\omega t} R_Q(t) dt \quad (18) $$

On voit ainsi que pour pouvoir accéder à des quantités telles que (15), (16), (17) ou (18), il faut savoir résoudre les problèmes (11) et (13), donnant $\rho$ et $p_S$. Pour les oscillateurs à 1 ddl décrits par des équations de type (2), la solution explicite du problème stationnaire (13) est connue dans des cas suffisamment généraux [10,11,12,19,20,22]. On dispose également de solutions explicites pour les oscillateurs non linéaires multidimensionnels, mais pour des situations beaucoup moins générales [19,20,22]. Par contre, on ne connaît pas à l'heure actuelle de solution analytique exacte pour le problème instationnaire (11) lorsque l'oscillateur est non linéaire, et ce même dans le cas scalaire (1 ddl).
Dans le cas général, il faut utiliser une méthode d'approximation pour résoudre (11) et le plus souvent (13) (notons que si l'on sait résoudre (11), la solution de (13) s'en déduit puisque, d'après les hypothèses, $p_S(x) = \lim_{t \to +\infty} \rho(x, t/x_0)$). Pour ce type d'équations aux dérivées partielles (parabolique pour (11), elliptique pour (13)), la méthode générale est la méthode de Galerkin à partir de laquelle on tire des procédures d'approximations variées en fonction du choix de la base fonctionnelle de projection. Deux principaux types de base fonctionnelle sont généralement utilisés. La première est construite à l'aide de polynômes orthogonaux avec poids tels, par exemple, que les polynômes d'Hermite [6,17,19,21]. La seconde consiste à utiliser la base des fonctions propres de l'opérateur $L^T$ ou $L$ (appelé générateur infinitésimal de la diffusion solution de (3)) [10,15,17,18].

Il faut toutefois noter que l'emploi de ces méthodes se heurte à deux difficultés importantes. D'une part, on est confronté aux problèmes généraux de vitesse de convergence et de troncature ; d'autre part, elles deviennent très difficiles d'utilisation, voire impraticables, lorsque la taille du problème (i.e. le nombre de ddl) augmente. En outre, ces difficultés sont encore accrues par le fait que l'opérateur de diffusion $a$ est toujours dégénéré (car associé à une EDSI obtenue à partir d'un système dynamique du second ordre).

Notons qu'il peut être très fructueux dans certains cas de remplacer le problème initial (3) par un problème approché du même type mais plus facile à résoudre. Un moyen pour y parvenir est par exemple d'appliquer la transformation de Van Der Pol à la diffusion solution de (3), puis d'utiliser la méthode de moyennisation stochastique [23,26,27,28,29,31]. On construit par ce procédé une approximation de la solution, qui est encore un processus de diffusion, mais dont la probabilité de transition et la loi invariante sont en général plus faciles à calculer.

## 4. LA MÉTHODE DES PERTURBATIONS ([23] à [31])

En toute rigueur, cette méthode ne s'applique que si l'équation dynamique non linéaire apparaît comme une perturbation faible d'une équation linéaire (ou, par extension, d'une équation non linéaire mais dont on connaît la solution). Notons qu'on peut souvent se ramener à cette situation par le biais de changements d'espace et de temps appropriés [23,31].

Pour fixer les idées, supposons que (2.1) puisse s'écrire sous la forme :

$$ \ddot{Q}(t) + 2\omega_0\xi_0\dot{Q}(t) + \omega_0^2 Q(t) + \varepsilon g(Q(t), \dot{Q}(t)) = D\dot{W}(t) \quad (19) $$

où $\omega_0$ et $\xi_0$ sont deux réels $> 0$, $g$ est une application non linéaire de $\mathbb{R}^2$ dans $\mathbb{R}$, dérivable autant de fois que nécessaire et vérifiant $g(-q, -\dot{q}) = -g(q, \dot{q})$, et $\varepsilon$ est un petit paramètre réel $> 0$. Pour simplifier, nous noterons encore $Q$ la solution stationnaire $Q_S$ de cette équation. Le principe de la méthode est alors le suivant. On suppose que $Q$ admet un développement asymptotique en puissances de $\varepsilon$, de la forme :

$$ Q(t) = Q_0(t) + \varepsilon Q_1(t) + \varepsilon^2 Q_2(t) + \dots \quad (20) $$
On introduit alors (20) dans (19) et on identifie les termes de même puissance en $\epsilon$. On obtient ainsi un système infini d'équations linéaires, qui s'écrit :

$$ \mathcal{L}_0 Q_0(t) = D \dot{W}(t) \qquad (21.1) $$

$$ \mathcal{L}_0 Q_1(t) = -g(Q_0(t), \dot{Q}_0(t)) \qquad (21.2) $$

$$ \mathcal{L}_0 Q_2(t) = -Q_1(t) \partial_q g(Q_0(t), \dot{Q}_0(t)) - \dot{Q}_1(t) \partial_q g(Q_0(t), \dot{Q}_0(t)) \qquad (21.3) $$

∴

où $\mathcal{L}_0$ est l'opérateur différentiel linéaire défini pour toute fonction $\varphi \in C^2(\mathbb{R}_+, \mathbb{R})$ par :

$$ \mathcal{L}_0 \varphi = \bar{\varphi} + 2\omega_0 \xi_0 \dot{\varphi} + \omega_0^2 \varphi $$

Les différentes équations de ce système se distinguent les unes des autres uniquement par leurs seconds membres, et le second membre de chacune d'elles n'est fonction que des solutions des équations qui la précédent. La résolution peut alors se faire en cascade. On remarquera que (21.1) correspond au problème linéaire associé à (19) ($\epsilon = 0$). Sa solution ($Q_0, \dot{Q}_0$) est donc bien connue : c'est un processus gaussien, stationnaire, centré, du second ordre, dérivable en moyenne quadratique, dont toutes les caractéristiques (fonction de corrélation, densité spectrale, lois finidimensionnelles) sont calculables explicitement (en résolvant l'EDSI linéaire associée à (21.1) ou plus directement en utilisant les techniques du filtrage linéaire des processus stationnaires).

Notant $h$ la fonction de réponse impulsionnelle causale associée à $\mathcal{L}_0$, les solutions stationnaires de (21) s'écrivent formellement :

$$ Q_0(t) = \int_0^\infty h(u) \dot{W}(t-u) du \qquad (22.1) $$

$$ Q_1(t) = - \int_0^\infty h(u) g(Q_0(t-u), \dot{Q}_0(t-u)) du \qquad (22.2) $$

∴

Introduisant (22) dans (20), on peut alors en déduire des approximations de quantités liées à la solution stationnaire $Q$ en tronquant le développement à un ordre désiré. On obtient, par exemple, au second ordre près, pour la fonction de corrélation $R_Q(\tau) = E[Q(t+\tau)Q(t)]$ et pour la densité spectrale $S_Q(\omega)$ :

$$ \bar{R}_Q(\tau) = R_{Q_0}(\tau) + \epsilon (R_{Q_0 Q_1}(\tau) + R_{Q_0 Q_1}(-\tau)) \qquad (23) $$

$$ \bar{S}_Q(\omega) = S_{Q_0}(\omega) + 2\epsilon \operatorname{Re}(S_{Q_0 Q_1}(\omega)) \qquad (24) $$

avec :

$R_{Q_0}(\tau) = E[Q_0(t+\tau)Q_0(t)]$ : fonction de corrélation de $Q_0$,

$R_{Q_0 Q_1}(\tau) = E[Q_0(t+\tau)Q_1(t)]$ : fonction d'intercorrélation de $Q_0$ et $Q_1$,

$S_{Q_0}(\omega) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} R_{Q_0}(\tau) e^{-i\omega\tau} d\tau$ : densité spectrale de $Q_0$,
$S_{Q_0Q_1}(\omega) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} R_{Q_0Q_1}(\tau)e^{-i\omega\tau}d\tau$ : densité interspectrale de $Q_0$ et $Q_1$, la fonction d'intercorrélation $R_{Q_0Q_1}(\tau)$ s'écrivant encore, en tenant compte de (22.1) et (22.2) :

$$R_{Q_0Q_1}(\tau) = \int_{-\infty}^{+\infty} E[Q_0(t+\tau)g(Q_0(t-u), \dot{Q}_0(t-u))]h(u)du \quad (25)$$

On voit ainsi que les différentes quantités intervenant dans (23) et (24) ne dépendent que du processus stationnaire $(Q_0, \dot{Q}_0)$ qui est entièrement déterminé. Bien entendu, on pourrait rechercher des approximations d'ordre plus élevé. Toutefois, pour les problèmes de grande dimension, l'expérience montre que, hormis le calcul du terme d'ordre 0 (correspondant au problème linéaire) et celui du terme d'ordre 1 (correspondant à la première contribution non linéaire), le calcul des termes d'ordre supérieur devient très vite impraticable. Cette méthode n'est donc véritablement intéressante que pour les faibles non linéarités. En outre, comme dans le cas déterministe, elle peut faire apparaître dans le développement asymptotique de la solution des termes séculaires préjudiciables à la qualité des approximations. C'est pourquoi on l'utilise généralement en lui associant la méthode de moyennisation stochastique [23,26,27,28,29,31] (et la transformation de Van Der Pol) qui a pour effet d'éviter l'apparition de tels termes. Signalons enfin que dans le cas stochastique, il n'existe pas, à notre connaissance, de résultats probants concernant la convergence de la série (20).

# 5. LES MÉTHODES DE LINÉARISATION ÉQUIVALENTE ([29] à [40])

## 5.1. La méthode classique

Cette méthode consiste à remplacer l'équation non linéaire (2.1) par une équation linéaire de la forme :

$$\ddot{Q}(t) + 2\omega_{eq}\xi_{eq}\dot{Q}(t) + \omega_{eq}^2 Q(t) = D\dot{W}(t) \quad (26)$$

qui soit «équivalente» à (2.1). Le mot «équivalent» a ici la signification suivante : réécrivons (2.1) sous la forme :

$$\ddot{Q}(t) + 2\omega_{eq}\xi_{eq}\dot{Q}(t) + \omega_{eq}^2 Q(t) + [f(Q(t), \dot{Q}(t)) - 2\omega_{eq}\xi_{eq}\dot{Q}(t) - \omega_{eq}^2 Q(t)] = D\dot{W}(t) \quad (27)$$

Nous dirons alors que (26) est stochastiquement équivalente à (2.1) si le terme entre crochets de (27), qui est un processus stationnaire, est minimum en moyenne quadratique : $E[(f(Q(t), \dot{Q}(t)) - 2\omega_{eq}\xi_{eq}\dot{Q}(t) - \omega_{eq}^2 Q(t))^2]$ minimum. Les coefficients $\omega_{eq}$ et $\xi_{eq}$ doivent donc vérifier les conditions du premier ordre :

$$\left\{ \begin{array}{l} \displaystyle \frac{\partial}{\partial \omega_{eq}} E[(f(Q(t), \dot{Q}(t)) - 2\omega_{eq} \xi_{eq} \dot{Q}(t) - \omega_{eq}^2 Q(t))^2] = 0 \\[1em] \displaystyle \frac{\partial}{\partial \xi_{eq}} E[(f(Q(t), \dot{Q}(t)) - 2\omega_{eq} \xi_{eq} \dot{Q}(t) - \omega_{eq}^2 Q(t))^2] = 0 \end{array} \right.$$
d'où l'on tire :

$$ \omega_{eq}^2 = E[Q(t)f(Q(t), \dot{Q}(t))] / E[Q^2(t)] \quad (28) $$

$$ 2\omega_{eq}\xi_{eq} = E[\dot{Q}(t)f(Q(t), \dot{Q}(t))] / E[\dot{Q}^2(t)] \quad (29) $$

Le processus figurant dans (28)-(29) est dans cette conception la solution stationnaire du système non linéaire. En toute rigueur, c'est donc par rapport à la loi de ce processus que doit être calculée l'espérance mathématique. Si l'on procède ainsi, on montre que la solution stationnaire du système linéarisé a exactement la même variance que celle du système non linéaire. On parle alors de linéarisation vraie ou en variance.

Toutefois, dans de nombreux cas, et notamment dans le cas multidimensionnel, on ne sait pas calculer la probabilité invariante du système non linéaire. On est alors amené à remplacer, dans les formules donnant les coefficients de la linéarisation, la loi invariante du système non linéaire par celle du système linéarisé, qui est gaussienne et que l'on sait entièrement caractériser. Dans ce cas, la linéarisation est qualifiée de gaussienne. On notera que cette deuxième approche nécessite la résolution d'un problème de point fixe car la loi invariante gaussienne du système linéarisé dépend des coefficients de la linéarisation, c'est-à-dire des inconnues du problème. En outre, elle ne fournit pas la variance exacte de la réponse stationnaire non linéaire, mais les approximations qu'elle en donne sont tout à fait correctes. Signalons enfin que l'hypothèse gaussienne permet, par une intégration par parties, et à condition que f soit dérivable, de donner à (28) et (29) la forme plus simple :

$$ \omega_{eq}^2 = E\left[\frac{\partial f}{\partial \dot{q}}(Q(t), \dot{Q}(t))\right] \quad (30) $$

$$ 2\omega_{eq}\xi_{eq} = E\left[\frac{\partial f}{\partial \dot{q}}(Q(t), \dot{Q}(t))\right] \quad (31) $$

D'une manière générale, la méthode de linéarisation classique (vraie ou gaussienne) présente les caractéristiques suivantes :

- elle reste constructive en grande dimension (notamment dans sa version «gaussienne»),

- elle donne des approximations de qualité pour les variances (voire les valeurs exactes dans sa version «vraie»),

- la solution linéarisée est toujours gaussienne (alors que la solution non linéaire ne l'est jamais) et sa densité spectrale est en général très éloignée de celle de la solution non linéaire.

En définitive, c'est une méthode qui présente l'avantage de rester opérationnelle en dimension supérieure, mais qui n'est vraiment intéressante que pour le calcul des variances. En fait, on l'utilise souvent comme étape intermédiaire dans d'autres méthodes d'approximation, et c'est sans doute dans cet emploi que réside son véritable intérêt.

## 5.2. Les méthodes de linéarisation à paramètres aléatoires [41,43,44]

Nous présentons ici deux méthodes exposées pour la première fois lors des journées sur le Contrôle Actif Vibro-Acoustique et la Dynamique Stochastique organisées par le L.M.A. du C.N.R.S. de Marseille les 9,10 et 11 octobre 1991 [41,44] et qui ont toutes deux pour vocation l'approximation des densités spectrales des solutions stationnaires d'équations dynamiques stochastiques non linéaires du second ordre.
La première est due à Ch.SOIZE de l'O.N.E.R.A. [44]. Celui-ci part d'une équation non linéaire du type (19), avec $\epsilon$ non nécessairement petit, dont il suppose a priori calculables les moments finis de tous ordres de la loi invariante de l'EDSI associée. Donc, si $X = (Q, \dot{Q})^t$ est la solution stationnaire de l'EDSI, il suppose calculables les moments $M_{\mu,\nu}$ de la v.a. $X(t) = (Q(t), \dot{Q}(t))^t$, tels que $M_{\mu,\nu} = E[Q(t)^\mu \dot{Q}(t)^\nu], |M_{\mu,\nu}| < +\infty$, où $\mu$ et $\nu$ sont des entiers $\ge 0$. Il cherche alors une équation linéaire équivalente à paramètres aléatoires de la forme :

$$ \ddot{Y}(t) + 2\omega_0 \xi_0 (1 + \epsilon \Theta) \dot{Y}(t) + \omega_0^2 (1 + \epsilon \Lambda) Y(t) = D \dot{W}(t) \quad (32) $$

où $(\Theta, \Lambda)$ est une v.a. à valeurs dans $\mathbb{R}^2$, indépendante de $W$ et de loi $P_{\Theta,\Lambda}(d\theta d\lambda)$ sur $(\mathbb{R}^2, B_2)$ devant vérifier certaines contraintes de support et d'intégrabilité.

Comme critère d'équivalence, SOIZE envisage d'abord le critère idéal, c'est-à-dire l'isonomie des solutions stationnaires de (19) et (32). Mais il montre très vite que le problème associé n'admet pas de solution. Il considère alors un critère d'égalité plus faible, portant uniquement sur les lois invariantes de (19) et (32). Mais là il se heurte à la difficulté de montrer l'existence et l'unicité de la solution du problème associé, qui se ramène dans ce cas à la résolution d'une équation intégrale. Cette question n'étant pas résolue à ce jour, il retient finalement un critère encore plus faible, mais viable celui-là, consistant en l'égalité des moments des lois invariantes, jusqu'à un certain ordre à choisir. Ainsi, si nous notons $M_{\mu,\nu}$ le moment d'ordre $\mu + \nu$ de la loi invariante de l'EDSI associée à (32), c'est-à-dire le moment $M_{\mu,\nu} = E[Y(t)^\mu \dot{Y}(t)^\nu]$ de la v.a. $Z(t) = (Y(t), \dot{Y}(t))^t$ où $Z = (Y, \dot{Y})^t$ est la solution stationnaire de l'EDSI associée à (32), ce critère s'écrit :

$$ M_{\mu,\nu} = M_{\mu,\nu}; (\mu,\nu) \in (0, \dots, \mu_{max}) \times (0, \dots, \nu_{max}) \quad (33) $$

où $\mu_{max}$ et $\nu_{max}$ sont deux entiers à choisir, tels que $|M_{\mu_{max},\nu_{max}}| < \infty$.

Pratiquement, cette méthode s'utilise comme suit. On se donne $\mu_{max}$ et $\nu_{max}$, ce qui définit la taille du système (33), et on calcule les $M_{\mu,\nu}$. On choisit ensuite une loi $P_{\Theta,\Lambda}(d\theta d\lambda)$ dans une famille de lois admissibles (i.e. vérifiant les contraintes de support et d'intégrabilité et comportant un nombre de paramètres compatibles avec la résolution du système (33)) et on calcule les $M_{\mu,\nu}$. Ce calcul se fait classiquement à partir de la loi conditionnelle, pour $\Theta$ et $\Lambda$ fixés, de la réponse stationnaire de (32) et de la loi $P_{\Theta,\Lambda}(d\theta d\lambda)$. Les moments $M_{\mu,\nu}$ dépendant des paramètres de cette loi, il ne reste plus alors qu'à résoudre le système non linéaire (33) pour déterminer ces paramètres. Notons qu'on pourrait également utiliser la démarche inverse en se donnant d'abord une loi admissible $P_{\Theta,\Lambda}(d\theta d\lambda)$ et en cherchant ensuite le nombre d'équations (33) à établir, c'est-à-dire le nombre de moments à retenir, pour pouvoir calculer les paramètres de la loi. Une fois la loi $P_{\Theta,\Lambda}(d\theta d\lambda)$ déterminée, le calcul de la densité spectrale $S_Y(\omega)$ de la solution stationnaire de (32) est immédiat. En effet, notant $S_{Y,\Theta}(\omega)$ la densité spectrale de la solution stationnaire de (32) pour $(\Theta, \Lambda)$ fixé à la valeur $(\theta, \lambda)$ (densité qui se calcule aisément car alors (32) est une équation linéaire à coefficients constants), $S_Y(\omega)$ s'obtient directement par :

$$ S_Y(\omega) = \int_{\mathbb{R} \times \mathbb{R}} S_{Y,\Theta,\lambda} P_{\Theta,\Lambda}(d\theta d\lambda) \quad (34) $$
Cette méthode a été testée par Ch.SOIZE sur plusieurs oscillateurs non linéaires de type Duffing [44] et les résultats qu'il obtient sont tout à fait satisfaisants. Nous l'avons nous-mêmes appliquée à un oscillateur de ce type décrit par l'équation :

$$ \ddot{Q}(t) + 2\omega_0 \xi_0 \dot{Q}(t) + \omega_0^2 Q(t) + \varepsilon Q(t)^3 = D \dot{W}(t) \quad (35) $$

avec les données numériques suivantes : $\xi_0 = 0,01$ ; $\omega_0 = 20\pi$ ; $\varepsilon = 78156566$, $D = 1$. La figure 1 montre les résultats obtenus en suivant exactement la démarche préconisée dans [44]. Ils font clairement apparaître que la méthode de linéarisation à paramètres aléatoires constitue un réel progrès par rapport à la méthode de linéarisation classique. Notons enfin que cette méthode semble devoir rester constructive en dimension supérieure.

<figure><img src="image_1.png" /><figcaption>Figure 1</figcaption></figure>

La seconde méthode est le fruit des travaux de R.BOUC du L.M.A. du C.N.R.S. de Marseille [41]. Pour présenter cette méthode, l'auteur part lui aussi d'une équation de type (19) (mais avec une non linéarité n'affectant que le déplacement $Q$, c'est-à-dire telle que $g(Q, \dot{Q}) = k(Q)$) qu'il commence par transformer à l'aide de changements d'espace et de temps appropriés, pour la ramener à la forme :

$$ \ddot{X}(t) + 2\xi_0 \dot{X}(t) + X(t) + \lambda f(X(t)) = \sqrt{4\xi_0} \dot{W}(t) \quad (36) $$

où le coefficient $\lambda$ ($> 0$) et la fonction $f$ sont déterminés par la transformation. Puis, s'appuyant sur des développements fondés sur l'emploi de la transformation de Van Der Pol (avec fréquence dépendant du processus enveloppe), de la méthode des perturbations et de
la méthode de moyennisation stochastique, il propose un modèle linéaire, stochastiquement équivalent à (36), de la forme :

$$ \ddot{Y}(t) + 2\xi_0 \dot{Y}(t) + \Omega^2(A) Y(t) = \sqrt{2\xi_0 [A\Omega(A)]^2} \dot{W}(t) \quad (37) $$

où A est une v.a. réelle indépendante de W, dont la loi $P_A(da)$ est la probabilité invariante d'une EDSI moyennisée, issue des développements mêmes de la méthode, et $\Omega^2$ est une fonction déterminée de $\mathbb{R}_+^*$ dans $\mathbb{R}$.

De ce modèle linéaire à paramètres aléatoires, il déduit alors facilement la densité spectrale de la réponse stationnaire, qui s'écrit :

$$ S_Y(\omega) = \int_{\mathbb{R}} S_{Y_a}(\omega) P_A(da) \quad (38) $$

où $S_{Y_a}$ est la densité spectrale de la solution stationnaire de (37) sous $A = a$ (densité qui se calcule sans difficulté puisque dans ce cas (37) est une équation linéaire à coefficients constants).

Tout comme Ch.SOIZE, R.BOUC a testé sa méthode sur des oscillateurs de Duffing et les approximations qu'il obtient des densités spectrales sont tout à fait remarquables (les deux premiers modes sont parfaitement décrits). On notera en outre, que rien ne semble s'opposer, du moins conceptuellement, à l'extension de cette méthode à la dimension supérieure.

## 6. CONCLUSION

Les méthodes que nous venons de présenter permettent toutes en principe d'obtenir des approximations des densités spectrales et des lois marginales des réponses stationnaires des systèmes dynamiques non linéaires du second ordre excités par des bruits blancs gaussiens.

Les méthodes de simulation donnent d'excellents résultats. Ce sont même les seules vraiment efficaces pour tous les problèmes bien posés. Mais l'importance des coûts numériques qu'elles engendrent limite pratiquement leur application aux problèmes de faible taille. On réserve donc généralement leur emploi à l'évaluation d'autres méthodes approxées. De plus, étant de nature purement numérique, on ne peut les utiliser pour étudier le comportement asymptotique des caractéristiques fonctionnelles estimées.

La méthode des équations de diffusion nécessite quant à elle la mise en ouvre de procédures d'approximation dont la construction soulève de nombreuses difficultés. Son emploi est très délicat et elle n'est effectivement utilisable qu'en faible dimension.

Le passage à la dimension supérieure constitue également un écueil pour la méthode des perturbations qui ne peut en pratique être utilisée, du moins dans sa conception classique, que pour les problèmes faiblement non linéaires.

La méthode de linéarisation stochastique classique présente deux avantages : d'une part, elle s'étend sans grandes difficultés à la dimension supérieure; d'autre part, elle fournit d'excellentes approximations des variances des solutions stationnaires. Mais, en contrepartie, elle ne permet pas d'obtenir de bonnes approximations des densités spectrales et des lois. En fait, cette méthode, comme d'ailleurs celle des perturbations, est souvent utilisée comme étape intermédiaire dans d'autres méthodes d'approximation.
Finalement, de toutes les méthodes présentées, seules les méthodes de linéarisation stochastique à paramètres aléatoires semblent en mesure de fournir une réponse au délicat problème de l'approximation explicite des densités spectrales et des lois. Bien entendu, elles n'en sont encore qu'à une phase préliminaire de leur développement, et elles demanderont sans doute à être améliorées, mais les résultats déjà obtenus par leurs auteurs en dimension un et pour des oscillateurs de Duffing sont très encourageants pour l'avenir.

Nous avons signalé en introduction deux raisons justifiant la nécessité de connaître les densités spectrales des réponses stationnaires des systèmes dynamiques stochastiques.

Bien entendu, ce ne sont pas les seules. Pour terminer, nous mettons l'accent sur l'une d'elles dont l'intérêt pratique nous paraît indéniable.

Avant cela, rappelons le résultat suivant (théorème de réalisation markovienne des processus gaussiens stationnaires physiquement réalisables [6,16]).

Soit $\omega \to S(\omega) : \mathbb{R} \to \text{Mat}_{\psi}(n, n)$ une densité spectrale et supposons satisfaites les hypothèses suivantes :

$$H1) \quad S \in L^1(\mathbb{R}; \text{Mat}_{\psi}(n, n))$$

$$H2) \quad S(\omega) \text{ est de rang plein pour tout } \omega$$

$$H3) \quad \int_{\mathbb{R}} \frac{\text{Log}(\det(S(\omega)))}{1 + \omega^2} d\omega > -\infty$$

$$H4) \quad \omega \to S(\omega) \text{ est rationnelle et admet la factorisation : } S(\omega) = \frac{R(i\omega) \cdot R^*(i\omega)}{|P(i\omega)|^2}$$

où $P(z)$ est un polynôme sur $\mathbb{C}$ à coefficients réels, de degré $q$, dont toutes les racines sont situées dans le demi-plan $Re(z) < 0$, et $R(z)$ est un polynôme sur $\mathbb{C}$ à coefficients dans $\text{Mat}_{\mathbb{R}}(n, n)$, de degré $r < q$.

Alors, sous ces hypothèses, il existe :

1) une matrice $A \in \text{Mat}_{\mathbb{R}}(m, m)$ asymptotiquement stable, avec $m = qn$,

2) une matrice $\Sigma \in \text{Mat}_{\mathbb{R}}(m, m)$, éventuellement dégénérée,

3) une matrice $B \in \text{Mat}_{\mathbb{R}}(n, m)$,

4) un processus de Wiener normalisé $W$ à valeurs dans $\mathbb{R}^m$,

5) une diffusion gaussienne homogène centrée $Y$ à valeurs dans $\mathbb{R}^m$,

6) un processus $Z$ gaussien, stationnaire, centré, à valeurs dans $\mathbb{R}^n$, admettant $S(\omega)$ pour densité spectrale,

tels que, à une isonomie près :

$$Z(t) = BY(t), \quad t \ge 0 \qquad (39)$$

avec $Y = (Y(t), t \in \mathbb{R}_+)$ solution stationnaire de l'EDSI linéaire sur $\mathbb{R}^m$ :

$$\begin{cases} dY(t) = AY(t)dt + \Sigma dW(t), & t > 0 \\ Y(0) = Y_0 \quad \text{p.s.} \end{cases} \qquad (40)$$

où $Y_0$ est une v.a. à valeurs dans $\mathbb{R}^m$, indépendante de $\{W(t), t \ge 0\}$ et ayant pour loi la probabilité invariante de (40-1).
Les matrices A, B et Σ se calculent directement à partir des termes de la factorisation de S(ω) en suivant une démarche standard issue de la démonstration même du théorème de réalisabilité [6,16]. L'entier m est appelé l'ordre de la markovianisation.

Revenons à présent au système général (1) et supposons, pour simplifier, que sa solution stationnaire $Q_S$ (dont nous admettons l'existence et l'unicité) soit centrée. Supposons en outre, hypothèse indispensable, que la densité spectrale $\omega \to S_{Q_S}(\omega) : \mathbb{R} \to \text{Mat}_{\oplus}(n, n)$ de cette solution soit connue.

Partant de là, observons que, si les hypothèses H1 à H3 peuvent être satisfaites par $S_{Q_S}(\omega)$, il n'en est pas de même de l'hypothèse H4. En effet, $S_{Q_S}(\omega)$ ne peut en aucun cas satisfaire cette dernière puisqu'elle ne peut déjà pas vérifier la propriété de rationnalité qui est une propriété propre aux densités spectrales des réponses stationnaires des systèmes linéaires.

Considérons alors une approximation $\tilde{S}_{Q_S}(\omega)$ de $S_{Q_S}(\omega)$ satisfaisant H4. Notons qu'il est toujours possible de construire une telle approximation. Le problème à résoudre, dit problème de factorisation spectrale fonctionnelle, s'énonce ainsi : trouver une fonction rationnelle matricielle de la forme $H(z) = R(z)/P(z)$, avec $R$ et $P$ satisfaisant H4, telle que $\tilde{S}_{Q_S}(\omega) = H(iw)H(iw)^*$ soit une approximation convenable (au sens de l'énergie) de $S_{Q_S}(\omega)$. La difficulté tient ici à la contrainte d'appartenance de $H$ à la classe de Hardy $H_+^2(\Pi^+, \text{Mat}_R(n, n))$, c'est-à-dire à la classe des fonctions analytiques sur $\Pi^+ = \{z/\text{Re}(z) > 0\}$ à valeurs dans $\text{Mat}_{\oplus}(n, n)$, qui sont les transformées de Laplace des fonctions de $L^2(\mathbb{R}_+, \text{Mat}_R(n, n))$. Une méthode de résolution est proposée dans [48,49] et donne d'excellents résultats dans le cas scalaire [1,51]. Citons également la méthode proposée dans [50] et appliquée dans [45,46], qui donne aussi de très bons résultats.

Supposons donc $\tilde{S}_{Q_S}(\omega)$ construite et supposons en outre qu'elle vérifie les hypothèses H1 à H3. Alors, d'après le théorème précédent, il existe un processus Z, gaussien, stationnaire, centré, défini par (39)-(40), dont la densité spectrale $S_Z$ coïncide avec $\tilde{S}_{Q_S}$ : $S_Z(\omega) = \tilde{S}_{Q_S}(\omega), \forall \omega \in \mathbb{R}$.

Donc, au sens de l'égalité des densités spectrales, et en ce sens seulement, le problème non linéaire (1) équivaut (à l'approximation au sens de l'énergie : $\tilde{S}_{Q_S}(\omega) \approx S_{Q_S}(\omega)$ près) au problème linéaire (39)-(40).

Ainsi, si l'on connaît la densité spectrale de la réponse stationnaire d'un système dynamique stochastique non linéaire, on peut construire un système dynamique stochastique linéaire à coefficients constants, spectralement équivalent au système non linéaire.

On imagine sans mal tout le parti à tirer d'une telle linéarisation. Elle permettrait notamment d'utiliser les codes de calcul spectraux existants (essentiellement valables pour les systèmes linéaires) pour traiter les systèmes non linéaires, une fois ceux-ci spectralement linéarisés.

D'où l'intérêt de disposer de méthodes permettant d'obtenir les densités spectrales (ou du moins d'en fournir de bonnes approximations) des réponses stationnaires des systèmes non linéaires.

Signalons que cette procédure de linéarisation, fait actuellement l'objet de développements dans les laboratoires de Génie Civil et de Mathématiques Appliquées de L'Université
Blaise Pascal de Clermont-Ferrand, dans le cadre d'une étude contractuelle financée par E.D.F.-S.E.P.T.E.N. [45,46,47].

## RÉFÉRENCES

[1] BERNARD P., FOGLI M., BRESSOLETTE Ph., LEMAIRE M. - Un algorithme de simulation stochastique par markovianasation approchée. Application à la mécanique aléatoire. *J. Méc. théor. et appl. (FRA)*, Vol.3, N.6, pp.905-950, 1984

[2] BERNARD P., FOGLI M., WAGNER S. - Un algorithme de simulation de processus gaussiens stationnaires par réalisation markovienne approchée. Application au calcul de la réponse d'une structure souple élancée soumise à une charge de vent. *Publications de l'A.F.R.E.M.*, Journées sur la Mécanique Aléatoire Appliquées à la Construction, pp. 157-165, juin 1984

[3] PRIESTLEY M.B. - Spectral analysis and time series. Academic Press Ed., New-York, 1981

[4] SHINOZUKA M. - Simulation of multivariate and multidimensional random processes. *J. of the Acoustical Society of America*, Vol.49, N.1, pp.357-367, 1971

[5] SHINOZUKA M., JAN C.M. - Digital simulation of random processes and its applications. *J. of Sound and Vibration*, Vol.25, N.1, pp.111-128, 1972

[6] SOIZE Ch. - Processus stochastiques et méthodes de résolution des problèmes aléatoires. *Cours de l'Ecole Centrale des Arts et Manufactures*, Paris, 1986

[7] SPANOS P.D. - Monte Carlo simulation of response of non-symmetric dynamic systems to random excitations. *Comps. Structs.*, Vol.13, pp.371-376, 1981

[8] TALAY D. - Simulation and numerical analysis of stochastic differential systems. *Publications de l'I.N.R.I.A. Sophia Antipolis*, 1989

[9] KETTANI-WAGNER S. - Etude comparative de techniques de simulations de fonctions aléatoires gaussiennes appliquées à l'analyse de la sécurité des structures. *Th. 3e cycle Génie Civil Clermont-Ferrand II*, mars 1986

[10] CAUGHEY T.K. - Non linear theory of random vibrations. *Advances in Applied Mechanic*, Vol.11, pp.209-253, 1971

[11] CAUGHEY T.K., MA F. - The exact steady state solution of a class of nonlinear stochastic systems. *Int. J. of Nonlinear Mechanics*, Vol.17, N.3, pp.137-142, 1983

[12] CAUGHEY T.K. - On the response of nonlinear oscillators to stochastic excitation. *Probabilistic Engineering Mechanics*, Vol.1, N.1, pp.2-4, 1986

[13] GUIKHMAN L., SKOROKHOD A. - The theory of stochastic processes. Springer Verlag Ed., Berlin, 1979

[14] IKEDA N., WATANABE S. - Stochastic differential equation and diffusion processes. *North Holland Ed.*, Amsterdam, 1981
[15] JOHNSON J.P., SCOTT R.A. - Extension of eigenfunction expansion solu-
tions of a Fokker-Planck equation. II. Second order system. *Int. J. of Nonlinear
Mechanics*, Vol.15, pp.41-56, 1980

[16] KREE P., SOIZE Ch. - Mécanique aléatoire. Dunod Ed., Paris, 1983

[17] RISKEN H. - The Fokker-Planck equation. Springer Verlag Ed., Berlin, 1989

[18] RYTER D. - On the eigenfunctions of the Fokker-Planck operator and of its adjoint. *Physica*, Vol.142A, pp.103-121, 1987

[19] SOIZE Ch. - Steady state solution of Fokker-Planck equation in higher dimension. *Probabilistic Engineering Mechanics*, Vol.3, N.4, pp.196-206, 1988

[20] SOIZE Ch. - Exact stationary response of multi-dimensional nonlinear hamiltonian dynamical systems under parametric and external stochastic excitations. *J. of Sound and Vibration*, Vol.149, N.1, pp.1-24, 1991

[21] WEN Y.K. - Approximate method for nonlinear random vibrations. *J. of the Engineering Mechanics Division*, EM4, pp.389-401, 1975

[22] ZHU W.Q., CAI G.Q., LIN Y.K. - On exact stationary response of stochas-
tically perturbed hamiltonian systems. *Probabilistic Engineering Mechanics*,
Vol.5, N.2, pp.84-87, 1990

[23] BOGOLIOUBOV N., MITROPOLSKI I. - Les méthodes asymptotiques en théorie des oscillations non linéaires. Gauthier-Villars Ed., Paris, 1962

[24] CRANDALL S.H. - Perturbation techniquess for random vibration of nonlinear systems. *J. of the Acoustical Society of America*, Vol.35, N.11, pp.1700-1705, 1963

[25] CRANDALL S.H., MARK D.W. - Random vibration mechanical systems. Academic Press Ed., New-York, 1973

[26] KHASMINSKII R.Z. - A limit theorem for the solutions of differential equations with random right-hand sides. *Theory of probability and its applications*, Vol.11, pp.390-405, 1966

[27] LIN Y.K. - Some observations on the stochastic averaging method. *Probabilistic Engineering Mechanics*, Vol.1, N.1, pp.23-27, 1986

[28] ROBERTS J.B., SPANOS P.D. - Stochastic averaging: an approximate method of solving random vibration problems. *Int. J. Nonlinear Mechanics*, Vol.21, N.2, pp.111-134, 1986

[29] ROBERTS J.B., SPANOS P.D. - Random vibration and statistical lineariza-
tion. John Wiley and Sons, New-York, 1990

[30] SOONG T.T. - Random differential equations in science and engineering. Academic Press, New-York, 1973
[31] **STRATONOVICH R.L.** – Topics in the theory of random noise. Gordon and Breach Ed., Vol.2, New-York, 1963

[32] **BERNARD P.** – About stochastic linearisation *IUTAM Symposium on Nonlinear Stochastic Mechanics*, Springer Verlag, 1991

[33] **BERNARD P., FOGLI M.** – Dynamique non linéaire des structures. Application à l'étude des vibrations aléatoires de structures tubulaires en présence de chocs. Rapport de contrat EDF SEPTEN/Université Blaise Pascal (Clermont-Ferrand), Vol.1, 1990

[34] **BRESSOLETTE Ph., FOGLI M.** – La méthode de linéarisation stochastique équivalente : application à l'étude des vibrations aléatoires d'un oscillateur scailaire à jeu. Rapport de recherche, Laboratoire de Génie Civil de l'Université Blaise Pascal de Clermont-Ferrand, sep.1991

[35] **CAUGHEY T.K.** – Equivalent linearization techniques. Journ. of Acoust. Soc. of Amer., Vol.35, N.11, pp.1706-1711, 1963

[36] **CRANDALL S.H.** – Heurisitic and equivalent linearization techniques for random vibrations of nonlinear oscillators. 6th. Int. Conf. Nonlinear Oscillations, Prague, 1978

[37] **IWAN W.D., MASON A.B.** – Equivalent linearization for systems subjected to non-stationary random excitation. Int. J. of Nonlinear Mechanics, Vol.15, pp.71-82, 1980

[38] **KOZIN F.** – The method of statistical linearization for nonlinear stochastic vibrations. *IUTAM Symposium, Innsbrück, Springer Verlag*, pp.45-56, 1988

[39] **SPANOS P.D.** – Stochastic linearization in structural dynamics. *Applied Mechanics Reviews*, ASME, Vol.34, N.1, pp.1-8, 1981

[40] **WEN Y.K.** – Equivalent linearization for hysteretic systems under random excitation. *J. of Applied Mechanics* ASME, Vol.47, pp.150-154, 1980

[41] **BOUC R.** – The power spectral density of a weakly damped strongly non-linear random oscillation, and stochastic averaging. *Publications du LMA-CNRS, Marseille, N.127*, pp.373-395, oct.1991

[42] **GUIHOT P.** – Analyse de la réponse de structures non linéaires sollicitées par des sources d'excitation aléatoires. Application au comportement des lignes de tuyauteries sous l'effet d'un séisme. *Thèse de Doctorat d'Université, spécialité Mécanique des Solides et des Structures, Université Paris VI*, 1990

[43] **MILES R.N.** – An approximate solution for the spectral response of Duffing's oscillator with random input. *J. of Sound and Vibration*, Vol.132, N.1, pp.43-49, 1989

[44] **SOIZE Ch.** – Sur le calcul des densités spectrales des réponses stationnaires pour des systèmes dynamiques stochastiques non linéaires. *Publications du LMA-CNRS, Marseille, N.127*, pp.297-344, oct.1991
[45] TAAZOUNT M. - Etude dynamique des structures à chocs sous sollicitations aléatoires. Thèse de Doctorat d'Université, spécialité Génie Civil, Université Blaise Pascal (Clermont-Ferrand), juin 1990

[46] BERNARD P., FOGLI M. - Etude des vibrations aléatoires d'une structure en présence de jeux : la méthode de linéarisation spectrale. Rapport de contrat EDF SEPTEN/Université Blaise Pascal (Clermont-Ferrand), Vol.2, 1991

[47] BERNARD P., BRESSOLETTE Ph., FOGLI M., TAAZOUNT M. - Oscillateurs à chocs : linéarisation spectrale et non linéarisation équivalente. Publications du LMA-CNRS, Marseille, N.127, pp.361-372, oct.1991

[48] BERNARD P., BONNEMOY C. - An algorithm for spectral factorization using random search techniques. Proba. Engng. Mech., Vol.4, N.2 pp.66-73, 1989

[49] BONNEMOY C. - Sur quelques aspects de l'utilisation de méthodes déterministes en milieu stochastique et inversement. Th. Doct. Sc. Math. Cl.Fd.II (FRA) sep-1987

[50] FRIOT E. - Simulation ARMA de processus stochastiques à partir d'une densité spectrale de puissance. Rapport de stage de fin d'études au LMA, Marseille, 1989

[51] BRESSOLETTE Ph. - Endommagement par fatigue de structures élancées en béton armé soumises au vent aléatoire. Th. Docteur-Ingénieur Génie Civil Clermont-Ferrand II, nov.1989

Laboratoire de Génie Civil
Université Blaise Pascal

63170 AUBIERE