Metoda Eulera-Maruyamy

Metoda Eulera-Maruyamy – metoda numerycznego rozwiązywania stochastycznych równań różniczkowych. Nazwa metody pochodzi od nazwiska japońskiego matematyka Gisiro Maruyamy (1916-1986), który rozszerzył metodę Eulera z 1768, stosowaną do numerycznego rozwiązywania równań różniczkowych zwyczajnych. Metoda ta jest jedną z niewielu metod numerycznych, które pozwalają metody numeryczne stosowane do równań deterministycznych rozszerzyć do równań stochastycznych^[1].

Opis metody

Równanie stochastyczne

Niech dane będzie stochastyczne równanie różniczkowe postaci

${\displaystyle \mathrm {d} X_{t}=a(X_{t})\,\mathrm {d} t+b(X_{t})\,\mathrm {d} W_{t}}$

oraz

${\displaystyle X_{T_{init}}=x_{0}}$ – zadany warunek początkowy (położenie cząstki w chwili początkowej ${\displaystyle T_{\mathrm {init} }}$),

gdzie:

${\displaystyle W_{t}}$ – proces Wienera,

${\displaystyle a(X_{t}),b(X_{t})}$ – zadane funkcje procesu losowego ${\displaystyle X_{t}.}$

Poszukiwane jest rozwiązanie tego równania na przedziale ${\textstyle [T_{init},T_{end}].}$

Analiza równania stochastycznego

O stochastycznym charakterze powyższego równania decyduje różniczka procesu losowego Wienera: gdyby była równa zeru, to proces byłby deterministyczny. Funkcja ${\displaystyle b,}$ stojąca przed tą różniczka, pełni zaś rolę wzmacniania lub osłabiania efektu fluktuacji losowych.

Z drugiej strony, funkcja ${\displaystyle a,}$ stojąca przed różniczką ${\displaystyle dt,}$ decyduje, na ile proces deterministyczny dominuje nad procesem losowym – gdyby funkcja ta była zerowa, to proces byłby czysto losowy.

Powyższe własności łatwo sprawdzić, rozwiązując podane dalej przykładowe równanie stochastyczne, zmieniając odpowiednio jego parametry.

Przybliżenie Eulera-Maruyamy

Przybliżeniem Eulera-Maruyamy dokładnego rozwiązania ${\displaystyle X_{t}}$ powyższego równania jest łańcuch Markowa ${\displaystyle Y_{n}}$ taki że:

• zdarzenia losowe łańcucha ${\displaystyle Y_{n},}$ ${\displaystyle n=0,1,2,\dots ,N,}$ zachodzą w dyskretnych chwilach czasu ${\displaystyle T_{\mathrm {init} }=\tau _{0}<\tau _{1}<\ldots <\tau _{N}=T_{\mathrm {end} },}$

odległych od siebie o wartość skokową

${\displaystyle \Delta t=(T_{\mathrm {end} }-T_{\mathrm {init} })/N,}$

• ${\displaystyle Y_{0}=x_{0}}$ – zadane zdarzenie początkowe,
• ${\displaystyle Y_{n+1},}$ ${\displaystyle n=0,1,2,\dots ,N-1}$ – zdarzenia obliczane wzorem rekurencyjnym:

${\displaystyle Y_{n+1}=Y_{n}+a(Y_{n})\,\Delta t+b(Y_{n})\,\Delta W_{n},}$

przy czym:

${\displaystyle \Delta W_{n}=W_{n+1}-W_{n},}$

${\displaystyle a(Y_{n}),b(Y_{n})}$ – zadane funkcje ${\displaystyle Y_{n}.}$

Dowodzi się, że zmienne losowe ${\displaystyle \Delta W_{n}}$ jako różnice zmiennych losowych o rozkładzie normalnym (proces Wienera) są niezależne oraz mają rozkład normalny o wartości oczekiwanej zero i wariancji ${\displaystyle \Delta t.}$

Przykład: Równanie Ornsteina-Uhlenbecka

Cztery losowe trajektorie procesu Ornsteina-Uhlenbecka dla θ = 1, σ = ${\displaystyle {\sqrt {2}}{:}}$
niebieska: położenie początkowe a = 10, μ = 0
pomarańczowa: położenie początkowe a = 0, μ = 0
zielona: położenie początkowe a = −10, μ = 0
czerwona: położenie początkowe a = 0, μ = −10
E(X) określa trajektorię, po której poruszałby się układ, gdyby czynniki losowe były zerowe (σ = 0).

Równanie stochastyczne

Niech dane będzie stochastyczne równanie różniczkowe jednej zmiennej ${\displaystyle X_{t}}$

${\displaystyle dX_{t}=\theta \cdot (\mu -X_{t})dt+\sigma \cdot dW_{t},}$

o warunku początkowym

${\displaystyle X_{0}=0,}$

gdzie ${\displaystyle \theta ,\mu ,\sigma }$ – stałe parametry; różniczka ${\displaystyle dW_{t}}$ odpowiada za proces losowy Wienera.

Powyższe równanie opisuje tzw. proces Ornsteina-Uhlenbecka.

Zamiana na postać dyskretna

Zgodnie z metodą Eulera-Maruyamy zamienia się to równanie na postać dyskretną

${\displaystyle \Delta Y_{n}=\theta \cdot (\mu -Y_{n})\cdot \,\Delta t+\sigma \cdot \Delta W_{n}.}$

Ponieważ ${\displaystyle \Delta Y_{n}=Y_{n+1}-Y_{n},}$ to otrzymuje się równania rekurencyjne

${\displaystyle Y_{n+1}=Y_{n}+\theta \cdot (\mu -Y_{n})\cdot \,\Delta t+\sigma \cdot \Delta W_{n},\qquad n=1,2,\dots ,N,}$

gdzie ${\displaystyle \Delta W_{n}}$ – zmienne losowe o rozkładzie normalnym i wariancji ${\displaystyle \Delta t.}$

Warunek początkowy przyjmie postać:

${\displaystyle Y_{0}=0.}$

Analiza zależność ruchu układu od parametrów równania

Zmieniając parametry równania stochastycznego, można otrzymać ruchu układu w zależności od różnych warunków otoczenia, z jakim układ oddziałuje. Np.

• przyjmując wartość parametru ${\displaystyle \sigma =0,}$ otrzyma się ruch czysto deterministyczny,
• przyjmując wartość parametru ${\displaystyle \theta =0,}$ otrzyma się ruch czysto losowy,
• parametr ${\displaystyle \mu }$ określa położenie, do którego zmierza układ po dłuższym czasie, jeżeli czynniki losowe nie są zbyt duże, tj. dla odpowiednio małej wartości ${\displaystyle \sigma$ ;} na rysunku pokazano cztery trajektorie układu w takich warunkach; mimo różnych wartości położeń początkowych układu (na rys. oznaczonych literą a), układ zmierza po pewnym czasie do położenia ${\displaystyle X=\mu .}$ Przykładem jest ruch wahadła poddanego działaniu losowego oddziaływania od otoczenia i jednocześnie tłumionego – niezależnie od położenia początkowego wahadło przyjmie ostatecznie najniższe położenie.

Kod programu

Poniżej podano przykład kodu w języku Python, całkujący równanie Ornsteina-Uhlenbecka. Program można testować, korzystając np. z darmowego notatnika colab google online.

Na ilustracji pokazano wynik symulacji komputerowej dla przyjętych w programie wartości parametrów.

Pięć przykładowych rozwiązań stochastycznego równania różniczkowego opisującego proces Ornsteina-Uhlenbecka, znalezionych numerycznie za pomocą metody Eulera-Maruyamy dla θ = 0.7, σ = 0.06, μ = 1.4 i identycznych położeń początkowych X= 0; na osi pionowej jest położenie układu, na osi poziomej czas. Każde rozwiązanie przedstawia możliwą, losowo przebiegającą trajektorią układu, która dla dłuższych czasów oscyluje wokół punktu X=μ =1.4

import numpy as np
import matplotlib.pyplot as plt

#Dane:
liczba_symulacji = 5
t_0 = 0    # chwila początkowa
t_1 = 7    # chwila końcowa
N   = 200  # liczba punktów podziału przedziału

theta  = 0.7
mu     = 1.4
sigma  = 0.06

#Definicja funkcji, która generuje liczbę losową przy każdym wywołaniu
def DW(Dt): 
    return np.random.normal(loc = 0.0, scale = np.sqrt(Dt))

# Część główna programu
Dt   = (t_1 - t_0) / N
t    = np.arange(t_0, t_1, Dt) # tablica dyskretnych chwil czasu
Y    = np.zeros(N)             # tablica dyskretnych położeń cząstki

for _ in range(liczba_symulacji):
    Y[0] = 0
    for i in range(0, N-1):
        Y[i+1]=Y[i] + theta*(mu -Y[i])*Dt +sigma*DW(Dt)
    plt.plot(t, Y)

plt.grid(True)  # Dodanie siatki do wykresu
plt.show()

Zobacz też

• równanie Langevina