Processus stochastiques et équations aux dérivés partielles

CHAPITRE IV :

Application

0. Introduction:

On applique dans ce chapitre quelques formules et méthodes citées dans les trois chapitres précédents. Cette application repose sur la simulation numérique. Pour commencer nous abordons la simulation des trajectoires de mouvement brownien qui est la base des autres simulations stochastiques de notre travail. Il faut discrétisé le problème, cette discrétisation se faisant au niveau des variables, par exemple sur la variable de temps pour le mouvement Brownien.

Après la simulation du mouvement Brownien, on procède à la simulation de processus de diffusion, que l'on utilise dans notre problème principal concernent l'interprétation probabiliste des EDPs.

1. Discrétisation du mouvement brownien :

Le mouvement Brownien standard est une variable aléatoire W(t) qui dépend continument de temps (t E [0, T]) et satisfait les quartes conditions connues(Chap1).

Pour la computation proposée il est utile de considérer le mouvement Brownien discrétisé, où W(t) est spécifie aux valeurs discrètes de t; ainsi nous donnons At = T /N pour un certain entier positive N et soit W~ une notation de W(ti) avec ti = jLt.

La première condition est W_o = 0 avec une probabilité 1, la 2^eme et la 3^eme sont données par :

Wi = Wi + AWi , j = 1,2, ...,N,

où OW~ est une variable aléatoire indépendante de -VAt.7V(0,1).

Le programme (PROG1) en MATLAB nous donne la simulation d'une trajectoire du Mouvement Brownien standard sur l'intervalle [0,1] avec N = 1000 ; le résultat de la simulation est illustré dans la figure (Fig.01) qui représente 1000 points (ti ,Bi) en les joignant par interpolation linéaire:

Remarque : tous les algorithmes des programmes (en MATLAB) cités dans ce chapitre sont donnés dans l'annexe B.

W(t)

-0.2

-0.4

0.8

0.6

0.4

0.2

1.4

1.2

0

1

simulation d'une trajectoire d'un mouvement brownien

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Fig 01: t

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5

Fig 02:

1

0.8

0.6

0.4

0.2

0

-0.2

-0.4

-0.6

-0.8

simulation d'une trajectoire d'un mouvement brownien dans 2D

Dans la figure (Fig.2) on donne la représentation d'une trajectoire de mouvement Brownien dans deux dimensions. On voit ici que le point de trajectoire marche aléatoirement dans tous les sens du plan.

Avec les mêmes démarches suivies dans la simulation d'une trajectoire d'un mouvement Brownien on peut simuler M trajectoires, le programme (PROG2) sur MATLAB, nous a donné la simulation de M = 1000 trajectoires, le résultat est illustré dans la figure(Fig.03) :

-2

-4

-5

-1

4

0

3

5

2

1

-3

W(t)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Fig.03: 1000 trajectoires de mouvement Brownien

t

2. Discrétisation d'un processus de diffusion:

2.1. Schéma d'Euler- Maruyama :

Une des simples méthodes de discrétisation de processus de diffusion est l'approximation d'Euler, appelée parfois Euler- Maruyama. On considère un

processus de diffusion X = (Xt, to t T} satisfait l'EDS :

dX_t = a( t, Xt)dt + b( t, X_t)dW_t (4.1)

dans [ to, T] avec la valeur initiale Xt_o = Xo.

Pour une discrétisation to = To
·
·
· TN = T de l'intervalle de temps [ to, T],

l'approximation d'Euler est un processus stochastique continu Y = (Y(t), to t T}
qui satisfait le schéma itératif :

Y_n+i = Y_n + c(r_n, Y_n)(r_n+i--T_n) + b(r_n, Y_n)(W_n+i--W_n) (4.2)

pour n = 0,1,2 ..., N -- 1 avec la valeur initiale Yo = X0, où nous avons écrit :

Yn = Y(rn)

Pour la valeur de l'approximation de temps de la discrétisation x_n, nous écrirons aussi : A_n= ;₁+1 -- x_n pour le n^tenie pas de temps et appelé S = max_n A_n le pas de temps maximale. Pour ce paragraphe on considère un temps de discrétisation équidistant x_n = to + nS avec S = A_nE A= (T -- to) /N pour N entier assez grand pour que DE [0,1] .

Quand le coefficient de diffusion est identiquement égale à zéro, l'itératif de schéma (4.2) réduit au schéma d'Euler déterministe pour l'équation différentielle ordinaire (EDO) :

dx

dt = a(t, x) (4.3)

La séquence (Y_n, n = 1,2., ..., N} des valeurs de l'approximation d'Euler (4.2) aux instants de temps discrétisé (x)8 = (r_n, n = 1,2, ..., N} peut avoir le même schéma

simulé dans le cas déterministe, la différence principale est que nous avons besoin à générer le pas aléatoire :

^AWn = Wrn-Ft--Wrn (4.4)

pour n = 0,1, 2 ..., N -- 1 de mouvement brownien W , on sait que ces pas sont des variables aléatoires gaussiens indépendants avec un moyenne E(AW_n) = 0 et de variance E(AW_n²) = A_n, donc on peut utiliser des nombres pseudo-aléatoires pour les incréments de mouvement brownien.

_n+i = Y_n + c/..A_n + b.AW_n

r

(4.5)

Y0 = X0

Pour simplifier nous écrivons : f = f(r_n,Y_n) pour une fonction définie sur IR⁺ x IR^l et n = 0, 1, 2 ..., N -- 1, alors nous pouvons écrire le schéma d'Euler (4.2) dans la forme abrévié :

La structure récurrente de schéma d'Euler qui évalue les valeurs approximatives du processus de diffusion aux instants de discrétisation est toujours la clé du Succès d'implémentation par ordinateur.

2.2. Exemple de discrétisation d'un processus de diffusion:

Nous illustrons l'aspect de simulation pour approximer le processus de diffusion, en examinant un simple exemple avec simple coefficient, on considère le processus X = (Xt, t 0) qui satisfait l'EDS linéaire :

dXt = aXtdt + bXtdWt (4.6)

Pour t E [0, 7] avec X0 E IR la valeur initiale, c'est un processus de diffusion avec drift a(t,x) = ax et un coefficient de diffusion b(t,x) = bx .

On peut avoir la solution exacte de cette EDS explicitement :

X_t = XoexpKa --¹1₂b²)t + bWt1 (4.7)

Pour t E [0, T] et W = (We, t 0) un Mouvement Brownien, connaître la solution de (4.6) explicitement donne la possibilité de comparer l'approximation d'Euler avec celle de la solution exacte, et calculer l'erreur.

Pour simuler la trajectoire de l'approximation d'Euler pour un temps de discrétisation, nous simplement donnons la valeur initiale Y0 = X0 et procéder par récurrence la prochaine valeur :

Y_n+i = Y_n + aY_nA_n + bY_nAW_n

pour n = 1,2, ... , accorder au schéma d'Euler avec drift et diffusion. Ici AW_n est la distribution gaussienne N(0, A_n) du mouvement brownien W sur le subintervalle x_n < t < T_n+1.

Pour la comparaison, nous pouvons utiliser (4.7) pour déterminer la valeur de la solution exacte pour un même échantillon de mouvement brownien, obtenu :

~

X,,, = Xoexp Ka -- 1 1₂ b²) r_n + b 1 Wi-il

~

Une application numérique de cet exemple en prenant l'intervalle de temps [0,1] , le pas A= 10^-2 et pour le processus de diffusion a = 0.5, b = 0.5 et X0 = 1. A travers le programme (PROG3) nous avons le résultat de la simulation de la solution exacte et la solution approximée par la méthode d'Euler de l'EDS (4.6). Le résultat est illustré dans la figure (Fig.04), où on peut voir que la trajectoire de solution approximée approche celle de la solution exacte.

solution exacte solution approximée

1.15

X(t)

1.1

1.05

1

0.95

0.9

0.85

0.8

0.75

0.7

0.65

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

t

Fig.04:solution exacte et approximée d'un processus de diffusion

3. Application à l'interprétation probabiliste des EDPs :

Dans cette partie on applique la méthode probabiliste sur quelques exemples des EDPs paraboliques, Pour voir l'efficacité de la méthode on applique en parallèle à ces EDPs une méthode déterministe (différences finies), on veut dire par efficacité, estce-que la méthode nous a donné une approximation acceptable de la solution ? On commence par poser le premier problème :

3.1. Problème N°=1 (équation de la chaleur):

On prend comme premier problème à traiter l'équation de la chaleur. Nous choisissons un problème que nous savons résoudre analytiquement et par différences finies pour essayer par une comparaison de voir l'efficacité et les avantages de la méthode probabiliste.

L'équation s'écrit comme suit :

lau(x, t) a ²u(x, 0

= a dans D = [0,T] x [0, xf]

at ax² (4.8)

u(x, 0) =

f(x)

avec :

1 -- domaine D = [0,0.1] x [0,1] et a = 1.

conditions aux limites

u(1,)

= 0

condition initiale f (x) = sinn-x . _fu(0, 0 = 0 (4.9)

On sait résoudre cette équation analytiquement et la solution (exacte) est :

u(x, t) = sin(n-x)exp(--n-²t)

La figure (Fig.05) présente la solution analytique u(x, t) en fonction de t et x:

Chapitre IV

Application

0.6

0.08

0.04

t

0

0

1

0.8

0.6

0.4

0.2

0

1

0.8

u(xt)

0.1

0.06

0.02

x

0.4

0.2

Fig.05:solution analytique en fonction de x et de t

3.1.1. Illustration numérique du problème par méthode déterministe (la méthode des différences finies) :

Cherchons l'approximation numérique de la solution u(x, t) qui satisfait l'équation (4.8) avec la condition initiale et les conditions aux limites (4.9). La méthode des différences finies nous permet de donner cette approximation ; on commence par discrétiser le problème. On divise le domaine de l'espace [0,xf] à M subintervalles, chacun de longueur Ax = xf/M, et on divise le domaine de temps [0,T] à N subintervalles, chacun de durée At = T/M, en utilisant la méthode d'Euler explicite (chap. II) on obtient :

k

ui₊1 - 214 +14_1

Ax²

k+1 k

ui -ui

= (4.10)

At

a

avecui^k

r-r^, u(xi,tk), la méthode explicite nous permet de résoudre le problème par

itération :

k+1 At

u_i = r(u+₁ + ut_₁) + (1 - 2r)u^k-, avec r = a

ex²

On a d'après le chapitre II que ce schéma est consistent, donc il suffit de choisir ?t et ?x pour que la condition de stabilité soit vérifiée, de cette façon on obtient un

schéma convergent. La condition de stabilité dans ce cas est que r ~ ~ .

~

Le programme (PROG.4) fournit une fonction qui calcule la solution à tous les points (tk,xt) explicitement, avec condition initiale et condition aux limites quelconques noté dans le programme par f(x) , bx0 et bxf :

Le programme (PROG.5) définit les fonctions : condition initiale et conditions aux limites (4.9), le pas de discrétisation ?t et ?x . Il suffit de donner par exemple M = 20, N = 300, donc r = 0.133 , et on trouve la solution:

La figure (Fig.06) donne l'illustration graphique de la solution approximée en fonction de t et x, par une première comparaison on voit la grande ressemblance entre la figure (Fig.05) et (Fig.06) :

On continue la comparaison par le calcul de l'erreur de troncature :

E = MaxtE~^k,1 < i < M,1 < k < N1

E = 1.4797 x 10^-4

0.06

0.04

0.2

0.4

x

t

0.02

0

0

Fig.06:solution approximée par la méthode des difference finies

1

0.8

u(x,t)

0.6

0.4

0.2

0

1

0.8

0.6

0.1

0.08

Puisque on connaît la solution exacte de l'équation (4.8), alors on peut calculer le vecteur erreur défini par :

I 1 k 1

et = k_.ei ,
·
·
· , _et ,
·
·
· , et

t )

avec :

e^k = lu~^k -- u(xi,tk)1, 1 < i < M, 1 < k < N

où u(xi, tk) est la solution analytique (exacte) au point (xi, tk) et ui^kest la solution approximée par la méthode des différences finies au même point.

La figure (Fig.07) illustre la solution exacte et approximée au long de l'axe t et pour un point d'espace fixé x = 0.5. On remarque que les deux graphes sont presque confondus ce qui veut signifier que l'erreur est acceptable; ce que montre la figure (Fig.08) qui illustre le vecteur d'erreur e(t, x) pour x = 0.5.

430

0.9

0.8

0.7

0.6

0.5

0.4

1

solution approximée solution exacte

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

t

Fig.07:solution exacte et approximée par la méthode des différences finies pour x=0.5

x 10^-4

e(t)

0.5

1.5

0

1

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

t

Fig.08: L'erreur de la méthode pour x=0.5

On calcule l'erreur quadratique pour xi = 1 / 2 par :

i1 VN

_e quad _ (e c)² (4.11)

_

N Li k=1 ^E

dans ce cas squad = 1.1699 x 10^-4

La Table.1 donne l'erreur maximum (troncature) et le temps (7'E) d'exécution de programme pour des résolutions différentes de l'équation (4.8) par la méthode des différences finies en changeant les valeurs de N et M (toujours pour x = 0.5):

N	90	600	1000	2000	8000	1.8x10^4
M	20	20	50	100	200	300
At	0.0013	1.66x10^-4	10^-4	5x10^-5	1.25x10^-5	5.5x10^-6
Ax	0.05	0.05	0.02	0.01	0.005	0.0033
r	0.44	0.066	0.25	0.5	0.5	0.5
ema x	0.0013	4.54x10^-4	6.05x10^-5	6.05x10^-5	1.51x10^-5	6.72x10^-6
7'E(s)	0.031	0.25	0.656	4.82	144.7	1076.53

Table.1

On constate de cette table que : chaque fois qu'on diminue la valeur de Ax il faut diminuer la valeur de At pour obtenir la stabilité de solution, cette situation provoque un système linéaire de grand dimension qui nécessite un temps d'exécution plus grand, mais l'erreur devient plus petite.

3.1.2. Illustration numérique du problème par la méthode probabiliste :

On passe par suite à l'application de la deuxième méthode sur le même problème, donc on s'intéresse à l'interprétation probabiliste de la solution de l'équation (4.8) pour t E [0, 7'], ce problème est traité dans le chapitre III, et par l'application de la formule de Feynman-Kac on peut écrire la solution sous forme :

u(x, t) = E _x[f(X_~)] (4.12)

où le processus stochastique sous-jacent (Xt) est la solution de :

X_t = x + f _o^t V2 dW_s ,0 t 7' (4.13)

La solution numérique est obtenue par la simulation de Monte-Carlo; l'application de schéma d'Euler donne une solution approximée de l'équation (4.13). On simule N trajectoires de cette solution, et par la fonction de condition initiale on calcule les valeurs f(X~), en moyennant ces valeurs on trouve :

~

u(x,t) ~ 1 ~

~~~

L'illustration graphique de cette solution est donnée par la figure (Fig.09) :

0

0

0.08

t

1.2

0.8

0.6

0.4

0.2

0

-0.2

1

x 0.8

0.6

0.4

0.2

0.1

0.06

0.04

0.02

1

u(x,t)

Fig.09: solution approximée par la méthode probabiliste

Le programme (PROG.6) donne la solution approximé par la méthode probabiliste, où le processus (X1 ~) est discrétisé avec même pas de temps utilisé à la première méthode (D.F), c.-à-d. on prend n = 300 et donc ?t = T /n = 1 / 300.

La figure (Fig.10) illustre la solution probabiliste et la solution exacte au long de l'axe de temps et pour un point d'espace fixé x = 0.5 , on peut voir que la solution probabiliste approche de la solution exacte avec des perturbations.

U(tpX)

0.9

0.8

0.7

0.6

0.5

0.4

1

solution probabiliste solution exacte

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

t

Fig.10:solution exacte et probabiliste pour x=0.5

On calcule l'erreur de troncature :

E = MaxtE ^k,1 k N)

E = 0.0227

L'erreur de troncature de cette méthode est plus grande que celle de la première. Ceci provient de l'ordre de convergence de la méthode de Monte Carlo.

Le vecteur d'erreur est donné par la formule :

e (xi, tk) = lUex(xi, tk) -- Uprob(xi, tk)I

où Uex est la solution exacte et Uprob est la solution par la méthode probabiliste.

La figure (Fig.11) illustre le vecteur e(x, t) pour une valeur de x fixé comme dans la première méthode x = 0.5.

e(tx)

0.025

0.015

0.005

0.02

0.01

0

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

Fig.11:l'erreur de la méthode probabiliste pour x=0.5 t

On calcule l'erreur quadratique : s qUad = 0.0098

On remarque que l'erreur quadratique est supérieure que celle de la première méthode. Ce qui signifie que la solution approximée par méthode des différences finies est meilleure que celle de la méthode probabiliste.

On peut affaiblir la valeur d'erreur par l'augmentation de la valeur de N (le nombre des trajectoires). La Table.2 donne l'erreur maximale des résultats de simulation de solution pour des valeurs différentes de N (toujours pour x = 0.5 ) :

N	500	1000	2000	8000	10000
emax	0.0110	0.0131	0.0098	0.0089	0.0065
TE (s)	0.141	0.734	3.313	215.21	330.82

Table.2

D'après la Table.2 on constate que la convergence de cette méthode est très lente. Mais le temps d'exécution de programme est inferieur à celle de la première méthode, car dans la méthode probabiliste on peut trouver la solution pour un x fixé directement, contrairement à la première méthode il faut trouver toute la solution (c'est à dire la matrice u(t,x)) et en suite extraire la solution pour un x fixé.

Interprétation des résultats :

Malgré que l'approximation de la solution par la méthode des différences finies soit meilleure que celle de la méthode probabiliste, on peut dire que la méthode probabiliste donne une approximation acceptable de la solution.

Puisque la méthode est basée sur la simulation de Monte Carlo, la convergence de la méthode est lente, et l'approximation sera faible de celles des méthodes déterministes.

L'utilisation de cette méthode est avantagée dans le cas de résolution des EDPs où la résolution par les méthodes déterministes implique une résolution des systèmes linéaires à grandes dimensions.

3.2. Problème N° 2 :

On prend comme deuxième problème l'équation aux dérivées partielles suivante:

{ au 1 0²u 1 au 1

x² (x, t)--x (x, t) + u(x, t) = 0 , dans D = [0,T] x IR,

at (x, t) + 2 ax² 2 ax ' 2 ' (4.14)

u(x, 0) = f(x)

{-- domaine D = [0,1] x [0,1] .

condition initiale f (x) = x² .

(4.15)

)u(0,t = 0

conditions aux limites {u(1, u(1, t) = exp (-- t / 2 )

On peut résoudre cette équation analytiquement et la solution (exacte) est :

u(x, t) = x²exp(--t/2)

La figure (Fig.12) illustre la solution en fonction de x et t :

Chapitre IV

Application

0.4

0.2

0.2

0.8

t

0

0

1

0.8

0.6

0.4

0.2

0

1

x

0.8

0.6

1

u(x,t)

0.4

0.6

Fig.12:solution analytique(exacte) en fonction de x et t

3.2.1. Illustration numérique de problème par la méthode des différences finies :

Cherchons l'approximation de la solution par la méthode des différences finies, en utilisant la méthode explicite (chap. II). On obtient le schéma suivant :

1

ur = 0-1B i -- r2Ai)u^li,__i + (1 -- 2 At -- riBi + r2Ai) ut -- r2Atu^li,₊1

avec r1 = At/Ax , r2 = AtiAx², Bt_. = -- ¹₂ x i , At = ¹₂ xi² .

On donne dans le programme (PROG.7) une fonction qui fait cette discrétisation et donne la solution approximée, la figure (Fig.13) représente cette solution en fonction de x et t.

Pour étudier l'erreur de la méthode, on fixe un point d'espace, comme dans le premier problème pour x = 0.5, la figure (Fig.14) illustre la solution exacte et la solution approximée par la méthode des différences finies, et la figure (Fig.15) illustre l'erreur commise par la méthode, d'après cette erreur on déduit que l'approximation est bonne, avec une erreur maximale emax = 1.5169 X 10~~

87

0

0

Chapitre IV

Application

0.6

0.4

0.2

0.4

0.2

0.8

t

0.6

x

1

0.8

0.6

0.4

0.2

0

1

0.8

u(x,t)

Fig.13:solution approximée par la méthode des différences finies

u(x,t)

0.19

0.18

0.17

0.16

0.15

0.14

0.13

0.12

0.11

0.2

0.1

solution approximée solution exacte

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 t 1

Fig.14:solution exacte et approximée par la méthode des différences finies pour x=0.5

1

e(x,t)

0.8

0.6

0.4

0.2

X 10^-4

1.6

1.4

1.2

0

1

Fig.15:le vecteur erreur pour x=0.5

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 t 1

3.2.2. Illustration numérique de problème par la méthode probabiliste :

La formule de Feynman-Kac permet de représenter la solution de problème :

1

u(x, t) = E_x [exp (- ₂ t) f(X_~)i , t E [0, T] (16)

avec f(x) = x² et le processus sous-jacent:

Xt = X0 -

1 i^t t 2J dt + f Xs dWs , X0 = x o o

Le programme (PROG.8) donne l'approximation de la solution en passant par la méthode de Monte Carlo :

La figure (Fig.16) représente la solution exacte et la solution approximée par la méthode probabiliste et la figure (Fig.17) représente l'erreur commise par cette méthode, avec N = 100 trajectoires.

L'erreur maximale emax = 0.0094, mais pour la plupart des points, l'erreur est inferieure à 0.005, avec une erreur absolue maximale égale à 5.9%, qui laisse une bonne approximation.

Interprétation des résultats :

Dans le deuxième problème l'importance de la méthode probabiliste est bien apparue, car l'application de la méthode des différences finies sur cet exemple est difficile à cause des coefficients non constant, c'est-à-dire, dépendent de la variable de l'espace, cette dépendance rend la solution non stable parfois, la méthode probabiliste grâce à la formule Feynman-Kac nous donne une bonne représentation de la solution, et plus on augmente le nombre des trajectoires dans la méthode de Monte Carlo on obtient une solution plus précise avec une erreur faible.

u(x,t)

0.26

0.24

0.22

0.18

0.16

0.2

solution exacte solution apprximée

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

t

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Fig.17: le vecteur erreur de la méthode

t

0.01

e(x,t)

0.009

0.008

0.007

0.006

0.005

0.004

0.003

0.002

0.001

0

Fig.16:solution exacte et solution approximée par la méthode probabiliste por x=0.5