Conception d'un pro logiciel interactif sous r pour la simulation de processus de diffusion

4.2.2 Modélisation d'une équation physique

Avant d'aborder la modélisation, nous annonçons deux théorèmes, donnant l'existe d'une relation simple entre les deux différentielles Itô et Stratonovitch, qui permet de passer de l'une à l'autre pour la résolution d'une équation différentielle stochastique

Théorème 4.1 (Différentielle d'Itô convertit en différentielle de Stratonovitch) Si un processus stochastique Xt satisfait l'équation d'Itô :

dXt = u(t,Xt)dt +a(t,Xt)dWt

alors il satisfait également l'équation de Stratonovitch :

dXt = u(t,Xt)dt + a(t,Xt)?dWt

oil le coefficient de dérive modifié u(t,Xt) est défini par :

u(t,x) = u(t,x) - 1 2 a(t,x)aa a(t,x)

x

Théorème 4.2 (Différentielle de Stratonovitch convertit en différentielle d'Itô) Si un processus stochastique Xt satisfait l'équation de Stratonovitch :

dXt = u(t,Xt)dt + ó(t,Xt) ?dWt

alors il satisfait également l'équation d'Itô :

dXt = u(t,Xt)dt +ó(t,Xt)dWt

oil le coefficient de dérive modifié u(t,Xt) est défini par:

u(t,x) = u(t,x) + ₂ó(t,x)^?ó(t,x)

1

?x

On considère une équation physique déterministe unidimensionnelle soumise à un bruit blanc î(t), sous la forme suivante

ÿx(t) = u(x,t)+ó(x,t)î(t) (4.8)

On suppose que le bruit est un processus centré E(î(t)) = 0 et de fonction de corrélation R(s,t) = E(î(s)î(t)) = ä(t - s). Le bruit blanc est modélisé par un processus de Wiener standard Wt. L'équation physique (4.8) se transcrit directement en équation de Stratonovitch

dXt = u(t,Xt)dt +ó(t,Xt)?dWt (4.9)

et se convertit en équation d'Itô sous la forme (appliquons le théorème 4.2)

dXt = A(t,Xt)dt +B(t,Xt)dWt (4.10)

Avec

A(x,t) = u(x,t) + ₂ó(x,t)^?ó(x,t)

1

?x

et

B(x,t) = ó(x,t)

D'où l'équation physique (4.8) se traduit en équation d'Itô sous la forme suivante

( )

dXt = u(t,Xt) + ¹ ₂ó(t,Xt)^?ó(t,Xt)dt+ó(t,Xt)dWt (4.11)

?x

L'équation de Fokker-Planck pour l'équation physique (4.8), s'écrit alors

?p ?t

t ~ ? u(x,t) · p + 1 ? = - 2ó(x,t) ? + 1 ?x (ó(x,t) · p) ?x2 (ó(x,t) · p) ?x 2

Dans le cas multidimensionnel, on note x = (x1,...,x_n). Le vecteur î(t) = (î1(t),...,î_m(t)) représente des bruits blancs centrés et de fonctions de corrélation normées

Ri j(s,t) = E(îi(t)îj(t)) = ä(t - s)

Les quantités u,ó1,...,ó_m sont des champs de vecteurs sur Rⁿ. L'équation physique multidimensionnel, s'écrit sous la forme

ÿx(t) = u(x,t) + ó1(x,t)î1(t) +
·
·
· + ó_m(x,t)î_m(t) (4.12)

se traduit mathématiquement par l'équation de Stratonovitch

dXt = u(t,Xt)dt +

m ? i=1

ói(t,Xt) 0 dWi(t)

oil W1(t),...,W_m(t) sont m processus de Wiener standard à valeurs dans Rⁿ. Cette équation correspond à l'équation d'Itô pour i = 1,...,n

1 ,---,^m dXl = u^i (t ,Xt) + , L z j=1 Si l'on pose

n ? k=1

j ?ai.(t,Xt)) m · (t x) ^i dt + ? ó^i(t ,Xt)dWj^i (t) , t ?xk j=1

n

?

i=1

L=

n

Ai(x) ? + 1 ? Bi j(x) ?2 ?xi 2 ?xi?xj

i, j=1

avec

?óⁱ _j(t,Xt)

ók j(t,Xt) ?xk

n

?

k=1

1 m

Aⁱ(x) = uⁱ(t,Xt)+ , ?

z j=1

et

n

B^ij(x) = ? óⁱk(x)ó^j_k(x)

k=1

alors le processus de diffusion Xt vérifie l'équation de Fokker-Planck

?p ?t

= Lp

Exemple 4.1 (L'oscillateur de Van Der Pol) On considère l'équation de Van Der Pol pour l'oscillateur ³ électrique, soumis à une fonce d'excitation aléatoire F, qui est centrée et de fonction de corrélation E(FtFt+h) = 2óä(h),

_X2

it + a (_b2 --11 ^Xÿ + ù²0X = F(t) (4.13)

L'équation de Van Der Pol (4.13) est utilisée pour modéliser des oscillateurs entretenus. Elle n'est pas linéaire et n'a pas de solution explicite. Les paramètres caractéristiques de cette équation sont : la pulsation ⁴ propre ù0, le paramètre de réaction ⁵ a et le paramètre de contrôle b.

3. Électricité : dispositif électrique qui sert à produire un courant alternatif périodique de fréquence déterminée.

4. Physique : augmentation momentanée et périodique de l'intensité d'une onde.

5. Physique : processus de modification de la structure d'un noyau atomique avec libération d'énergie.

On écrit cette équation du second degré sous la forme d'un système d'équations du premier degré

(X^ÿ = Y(X2 ) (4.14)
Yÿ = -a b2 - 1 Y - ù² ₀X + F(t) se traduit mathématiquement par l'équation de Stratonovitch

(

dXt = Ytdt ( _(X2 ) )

dYt = - a b2 - 1 Yt + ù2 0Xt dt + 2ó ? dWt

l'équation (4.13) se traduit à un système d'équations d'Itô sous la forme suivante

IdXt = Ytdt( (X2 ) ) (4.15)
dYt = - a b2 - 1 Yt + ù2 0Xt dt + 2ódWt L'équation de Fokker-Planck pour l'équation de Van Der Pol (4.13), s'écrit alors

?p
?t

~ _'x2 ~ ~

? ? + ? 0xp) + 1 ?2

= - _?x(yp) + a b2 - 1 yp ?y^(ù2 ?y2 (2óp)

?y 2

= -y	~x2 ~ ? ? ?y p + ó ?2 ?x p + a b2 - 1 ?y(yp) + ù^2 0x ? ?y2 p

Simulation numérique de l'équation de Van Der Pol

La fonction snssde2D⁶ permettre de simulée numériquement la solution approchée des systèmes d'équations différentielles stochastiques deux dimensions par des schémas numériques classiques.

R> help("snssde2D")

R> example("snssde2D")

R> snssde2D(N, T = 1, t0, x0, y0, Dt, driftx, drifty, diffx, diffy,

+ Step = FALSE, Output = FALSE, Methods = c("SchEuler",

+ "SchMilstein", "SchMilsteinS", "SchTaylor", "SchHeun",

+ "SchRK3"), ...)

6. simulation numerical solution of stochastic differential equations two-dimensional.

Details:

N La taille de processus.

M Le nombre de trajectoire à simulée.

T L'instant final (par défaut égale à 1).

t0 L'instant initial.

x0 & y0 Les valeurs initiaux.

Dt La discrétisation où le pas (Si Dt est fixée alors par défaut

T = t0 + Dt * N).

driftx & drifty Coefficient de dérive (une expression qui dépende de x, y et de t).

diffx & diffy Coefficient de diffusion (une expression qui dépende de x, y et de t).

Output Pour sauvegardée les résultats de simulation sous forme Excel

(par défaut c'est FALSE).

Step Pour l'animation, simulation étape par étape (par défaut c'est FALSE).

Methods Différents méthodes de simulation (par défaut schéma d'Euler),

avec SchEuler(3.23), SchMilstein(3.24),SchMilsteinS(3.25), SchTaylor(3.26), SchHeun(3.27), SchRK3(3.28).

Posant par exemple (b,ù0,ó) = (4,2,0.5) et x0 = y0 = 0. On remarque deux états du ce système physique très différents pour le paramètre a = 0 et a > 0,donc on a l'équation physique

f X· +4X = F(t), a = 0

~X2 ~

X· + 2 16 _ 1 ^Xÿ + 4X = F(t), a > 0

se convertit mathématique en équation d'Itô

(

dXt = Ytdt , a = 0

dYt = _4Xtdt + dWt

et

(

dXt = Ytdt ( _(X2 ) ) , a > 0

dYt = _ 2 16 _ 1 Yt + 4Xt dt + dWt

R> fx <- expression( y ) R> gx <- expression( 0 )

R> fy1<- expression( -4*x ) R> gy <- expression( 1 )

R> snssde2D(N = 10000, T = 1, t0 = 0, x0 = 0, y0 = 0, Dt = 0.1,

+ driftx = fx, drifty = fy1, diffx = gx, diffy = gy)

R> fy2<- expression( -(2*((x/4)^{^}2 -1)*y + 4*x) )

R> snssde2D(N = 10000, T = 1, t0 = 0, x0 = 0, y0 = 0, Dt = 0.1,

+ driftx = fx, drifty = fy2, diffx = gx, diffy = gy)

FIGURE 4.1 - L'oscillateur de Van Der Pol, régime permanent sinusoïdal a = 0.

FIGURE 4.2 - L'oscillateur de Van Der Pol, régime permanent non sinusoïdal a > 0.