Projet : OMEGA

Précédent : Fondements scientifiques Remonter : Fondements scientifiques Suivant : Analyse stochastique des algorithmes

   
Méthodes probabilistes pour les équations aux dérivées partielles

Participants : Mireille Bossy, Madalina Deaconu, Bernard Roynette, Denis Talay, Pierre Vallois.

De nombreux problèmes d'évolution linéaires ou non linéaires

 

$${\frac{du}{dt}}$$ (1)

peuvent être interprétés à l'aide de processus de Markov bien choisis : on interprète u à l'aide du générateur infinitésimal du semigroupe de transition d'un processus de Markov (X_t) ou bien à l'aide de l'adjoint de ce générateur. Les motivations de cette démarche peuvent être d'ordre théorique et/ou numérique. En effet, en particulier lorsque X = (X_t) est solution d'une équation différentielle stochastique, le calcul stochastique permet parfois d'obtenir des résultats d'existence, d'unicité ou de régularité de la solution de (1) plus efficacement que les techniques d'analyse habituelles : le théorème de Girsanov, le calcul de Malliavin, la propagation du chaos sont des outils puissants qui n'ont pas d'analogues en analyse « déterministe » des équations aux dérivées partielles. D'autre part, dès que l'on peut écrire la solution de (1) sous la forme d'une espérance du type u(t) = EF(X ^_. ) avec F fonctionnelle sur l'espace des trajectoires de X entre 0 et t, on peut chercher à développer une méthode de Monte-Carlo pour approcher u(t) même si on ne sait pas simuler des trajectoires exactes de X : il suffit de construire un processus proche (en loi) de X, en simuler un grand nombre de trajectoires entre 0 et t, évaluer la fonctionnelle F le long de chaque trajectoire simulée et enfin moyenner toutes les valeurs obtenues.

Donnons un exemple élémentaire. Considérons l'équation de la chaleur

 

$$\frac{\partial u}{\partial t}(t,x)=\nu\Delta u(t,x),~\forall (t,x)\in]0,T]\times\mathbbm {R}^d$$ (2)

avec pour condition initiale u(0, ^. ) = u₀( ^. ) une fonction mesurable bornée. Le paramètre $\nu$ est strictement positif, et est appelé « paramètre de viscosité » en mécanique des fluides ou « volatilité » en finance.

On vérifie facilement que la fonction $$\forall(t,x)\in]0,T]\times\mathbbm {R}^d,~u(t,x):=\mathbbm {E}u_0(x+\sqrt{2\nu} W_t)$$où (W_t) est un mouvement brownien standard à valeurs dans $\mathbbm {R}^d$ satisfait (2) ainsi que $\lim_{t\rightarrow0}^{}$u(t, x) = u₀(x) en tout point de continuité de u₀. Par application de la loi des grands nombres, on peut donc approcher u(t, x) par

$${\frac{1}{N}}$$$$\sum_{i=1}^{N}$$u₀(x + $$\sqrt{2\nu t}$$g_i($$\omega$$))

où les {g_i($\omega$)} forment une famille de variables aléatoires gaussiennes indépendantes, à valeurs dans $\mathbbm {R}^d$, centrées et de matrice de covariance $Id_{\mathbbm {R}^d}$. Cet algorithme est très simple à mettre en oeuvre : on sait effectuer des tirages gaussiens indépendants à l'aide d'appels à un générateur de nombres pseudo-aléatoires uniformément répartis ; en outre il est naturellement parallélisable : le i^e processeur a la tâche d'engendrer g_i($\omega$). La vitesse de convergence est décrite par des théorèmes-limite tels que le théorème de limite centrale, la loi du logarithme itéré, l'inégalité de Berry-Esseen : la convergence est d'ordre 1/$\sqrt{N}$, elle est donc lente. Toutefois, le coût de l'algorithme croît seulement linéairement avec la dimension d de l'espace puisqu'on simule Nd trajectoires d'un mouvement brownien unidimensionnel standard, et ce coût est indépendant du paramètre $\nu$.

Typiquement, les méthodes de Monte-Carlo pour des équations aux dérivées partielles elliptiques ou paraboliques peuvent permettre de traiter des problèmes extrêmes, en très grande dimension ou avec de très faibles viscosités, lorsqu'il serait difficile, ou démesurément coûteux, d'utiliser des algorithmes classiques.

Soit à présent le problème parabolique

 

$${\frac{\partial u}{\partial t}}$$(t, x) = $$\sum_{i=1}^{d}$$bi(t, x)$${\frac{\partial u}{\partial x_i}}$$(t, x) + $${\textstyle\frac{1}{2}}$$$$\sum_{i,j=1}^{d}$$aij(t, x)$${\frac{\partial^2 u}{\partial x_i\partial x_j}}$$(t, x),
$$\forall$$(t, x) $$\in$$ ]0, T] x Rd,

(3)

où b est une fonction à valeurs dans $\mathbbm {R}^d$ et a une fonction à valeurs dans l'espace des matrices symétriques définies positives. Sous certaines conditions, on sait que l'unique solution régulière vérifie $$u(t,x)=\mathbbm {E}u_0(X_T^{t,x}),$$où (X_$\theta$^{t, x}) est le processus de Markov solution d'une équation différentielle stochastique

 

$$\theta \int_{t}^{\theta} \sum_{j=1}^{r} \int_{t}^{\theta} \sigma_{j}^{}$$ (4)

où les matrices $\sigma$(t, x) sont des racines carrées des matrices a(t, x). La discrétisation en temps de (4) conduit naturellement au processus de Markov à temps discret suivant :

 

$$\left\{\vphantom{ \begin{array}{lcl} \bar X^{t,x}_0 &=& x, \ \b... ...{n}}^{t,x})(W^j_{(p+1)\frac{T-t}{n}}-W^j_{p\frac{T-t}{n}}). \end{array}}\right. \begin{array}{lcl} \bar X^{t,x}_0 &=& x, \ \bar X_{(p+1)\frac{T... ...rac{T-t}{n}}^{t,x})(W^j_{(p+1)\frac{T-t}{n}}-W^j_{p\frac{T-t}{n}}). \end{array} \left.\vphantom{ \begin{array}{lcl} \bar X^{t,x}_0 &=& x, \ \ba... ...{n}}^{t,x})(W^j_{(p+1)\frac{T-t}{n}}-W^j_{p\frac{T-t}{n}}). \end{array}}\right.$$ (5)

Il est facile de simuler des trajectoires indépendantes ($\bar{X}_{\cdot}^{t,x}$($\omega_{i}^{}$)) de ce processus puisque les variables aléatoires W^j_{(p + 1)(T - t)/n} - W^j_{p(T - t)/n} sont mutuellement indépendantes et de même loi gaussienne centrée de variance (T - t)/n. On peut donc numériquement approcher u(t, x) par

u(t, x) $$\sim$$ $${\frac{1}{N}}$$u₀($$\bar{X}_{T}^{t,x}$$($$\omega_{i}^{}$$)).

La vitesse de convergence de la méthode dépend à la fois du nombre Nde simulations et du nombre n de pas de temps.

Le procédé s'étend dans des directions variées : problèmes elliptiques, problèmes de transport (applications en neutronique), problèmes avec conditions frontière de Dirichlet ou de Neumann, problèmes intégro-différentiels, etc.

Au lieu de vouloir résoudre (3), on peut s'intéresser au problème adjoint

 

$${\frac{\partial p}{\partial t}}$$(t, x) = - $$\sum_{i=1}^{d}$$$${\frac{\partial}{\partial x_i}}$$(bi(t, x)p(t, x)) + $${\textstyle\frac{1}{2}}$$$$\sum_{i,j=1}^{d}$$$${\frac{\partial^2}{\partial x_i\partial x_j}}$$(aij(t, x)p(t, x)),
$$\forall$$(t, x) $$\in$$ ]0, T] x Rd,

(6)

avec la condition : p(t, x)dx converge faiblement vers une mesure de probabilité donnée lorsque t tend vers 0. Supposons (ce n'est pas une restriction) que cette mesure soit la masse de Dirac en x. Soit (X ^{_. 0, x}($\omega_{i}^{}$)) des trajectoires indépendantes de la solution de (4) avec t = 0. Sous de bonnes hypothèses, la mesure empirique

$$\mu_{N}^{}$$ : = $${\frac{1}{N}}$$$$\sum_{i=1}^{N}$$$$\delta_{X_t^{0,x}}^{}$$($$\omega_{i}^{}$$)

converge faiblement vers p(t, x)dx. Cette remarque sous-tend une famille de méthodes particulaires stochastiques pour les équations aux dérivées partielles non linéaires de type équation de McKean-Vlasov :

 

$$\left\{\vphantom{ \begin{array}{ll} \displaystyle{\frac{\partial ... ...,~(x,t) \in\mathbbm {R}^d\times]0,T],\ U_{t=0}\:=\:U_0. \end{array} }\right. \begin{array}{ll} \displaystyle{\frac{\partial U_t}{\partial t}}=... ...y)\right),~(x,t) \in\mathbbm {R}^d\times]0,T],\ U_{t=0}\:=\:U_0. \end{array} \left.\vphantom{ \begin{array}{ll} \displaystyle{\frac{\partial U... ...,~(x,t) \in\mathbbm {R}^d\times]0,T],\ U_{t=0}\:=\:U_0. \end{array} }\right.$$ (7)

La fonction b( ^. , ^. ) à valeurs dans $\mathbbm {R}^d$ qui intervient dans la partie non linéaire de l'équation est appelée noyau d'interaction. L'équation ci-dessus est considérée au sens des distributions. La théorie probabiliste de la propagation du chaos montre que la solution U_t s'interprète à l'aide de la loi limite d'un système de particules interagissant entre elles. La dynamique des particules est décrite par le système différentiel stochastique de dimension N x d

 

$$\left\{\vphantom{ \begin{array}{ll} X^{i}_t =\displaystyle{\int_0... ...ire de loi } U_0, ~\hbox{indépendante de } X^j_0,~i\not=j. \end{array} }\right. \begin{array}{ll} X^{i}_t =\displaystyle{\int_0^t}\displaystyle{\... ...le aléatoire de loi } U_0, ~\hbox{indépendante de } X^j_0,~i\not=j. \end{array} \left.\vphantom{ \begin{array}{ll} X^{i}_t =\displaystyle{\int_0^... ...ire de loi } U_0, ~\hbox{indépendante de } X^j_0,~i\not=j. \end{array} }\right.$$ (8)

La propagation du chaos implique la convergence au sens des mesures, quand N tend vers l'infini, de la mesure empirique 1/N$\sum_{i=1}^{N}$$\delta_{X^{i}_t}^{}$ vers U_t. En particulier, un lissage par convolution de la mesure empirique converge vers la fonction U_t. À partir de cette interprétation probabiliste, on développe un algorithme d'approximation de U_t fondé sur la simulation du système de particules (Xⁱ_t, 1 $\leq$ i $\leq$ N) ; la mesure initiale U₀ est approchée par une combinaison linéaire de masses de Dirac, ce qui fournit les positions initiales des particules, qu'on déplace en simulant une (et une seule) réalisation approchée du système (Xⁱ_t, 1 $\leq$ i $\leq$ N) ci-dessus.

La complexité de l'analyse de la vitesse de convergence dépend essentiellement de la singularité éventuelle du noyau d'interaction b( ^. , ^. ). Pour la plupart des équations provenant de problèmes physiques (et en particulier pour les équations de Burgers ou de Navier-Stokes en dimension 2), le noyau d'interaction est singulier. La vitesse de convergence dépend du nombre N de particules et du pas de temps utilisé pour la discrétisation de (8).

Pour un aperçu de résultats sur les méthodes de Monte-Carlo et certaines méthodes particulaires stochastiques, on pourra consulter [[8]].

Footnotes

...$\mathbbm {R}^d$

Un processus stochastique est une famille de variables aléatoires sur un espace probabilisé $(\Omega,{\cal F},\mathbbm {P})$indicées par le temps : $\{X_t(\omega),t\in\mathbbm {R}^+\}$ ; à $\omega$ fixé l'application t $\rightarrow$ X_t($\omega$) est appelée « trajectoire ». Un exemple de processus est le mouvement brownien standard à valeurs dans $\mathbbm {R}^d$, processus à trajectoires presque sûrement continues défini de la manière suivante : W₀ = 0 presque sûrement ; pour tout 0 < s < t, la variable aléatoire W_t - W_s est de loi gaussienne centrée, de matrice de covariance $(t-s)Id_{\mathbbm {R}^d}$, indépendante de la famille {W_$\theta$, 0 $\leq$ $\theta$ $\leq$ s}.

Précédent : Fondements scientifiques Remonter : Fondements scientifiques Suivant : Analyse stochastique des algorithmes