Z Skrypty dla studentów Ekonofizyki UPGOW
(Różnice między wersjami)
|
|
Linia 61: |
Linia 61: |
| | | |
| <math> X^j(t_i) = X^j(t_{i-1}) + F^j (\mathbf X(t_{i-1}), t_{i-1}) h + \sqrt{h} \sum_{k=1}^{n} G^{j,k}(\mathbf X(t), t) N^k(0,1)\; </math> | | <math> X^j(t_i) = X^j(t_{i-1}) + F^j (\mathbf X(t_{i-1}), t_{i-1}) h + \sqrt{h} \sum_{k=1}^{n} G^{j,k}(\mathbf X(t), t) N^k(0,1)\; </math> |
- |
| |
- | === Proces Wienera ===
| |
- |
| |
- | Proces Wienera jest rozwiązaniem następującego stochastycznego równania różniczkowego:
| |
- |
| |
- | <math>dX(t)= \sqrt{2 D} dW(t)\;</math>.
| |
- |
| |
- |
| |
- | Jego realizacja jest funkcją ciągłą, ale jednocześnie nigdzie nie jest różniczkowalna.
| |
- | [[Plik:wiener.png|thumb|360px|Proces Wienera]]
| |
- | Stosując schemat Eulera można wygenerować pojedyńczą trajektorię takiego procesu:
| |
- | <source lang='matlab'>
| |
- | h=0.01;
| |
- | N=400;
| |
- | x(1)=0;
| |
- | D=1;
| |
- | for i=2:N
| |
- | x(i)=x(i-1)+ sqrt(2*D*h)* normrnd (0,1,1,1);
| |
- | endfor
| |
- |
| |
- | plot((1:N)*h,x,'-')
| |
- | </source>
| |
- |
| |
- | Jest to trajektoria, która w granicy <math>h\to\infty</math> jest funkcją nigdzie nie różniczkowalną. Numerycznie przejawia się to w tym, że niezależnie od tego jak małe h weźmiemy do symulacji, wykres procesu Wienera zawsze będzie wyglądał na "zaszumiony". Można też korzystając z techniki Mostu Browna z danego przybliżenia procesu Wienera dla pewnego h wysymulowac przybliżenie dla h/2. Również wtedy można stwierdzić, że wykres nigdy nie będzie wizualnie "gładki", co w tym przypadku oznacza matematycznie nieciągłość w każdym punkcie.
| |
- |
| |
- | === Proces Wienera: rozkład P(x,t) ===
| |
- | Proces Wienera (symetryczny) spełnia, jak wiadomo równanie dyfuzji:
| |
- |
| |
- | <math>\frac{\partial P(x, t)}{\partial t} = D \frac{\partial^2 P(x, t)}{\partial x^2} </math>
| |
- |
| |
- | Rozważmy przypadek w którym cząstka staruje w punkcie <math>x=0</math> w czasie <math>t=0</math>. W języku gęstości prawdopodobieństwa oznacza to
| |
- |
| |
- | <math>P(x, 0) = \delta(x)\;</math>
| |
- |
| |
- | Jego rozwiązaniem na równania dyfuzji prostej z takim warunkiem początkowym jest rozkład Gaussa:
| |
- |
| |
- | <math>P(x, t) = \frac{1}{\sqrt{4\pi D t} }\; \exp \left[ - \frac{x^2}{4Dt}\right]</math>
| |
- |
| |
- | w którym odchylenie standardowe jest proporcjonalne do czasu.
| |
- | Zweryfikujmy ten fakt numerycznie. W tym celu wysymulujemy trajektorie 20000 procesów Wienera, dla 100 kroków. Poprzedni program łatwo można zmodyfikować by obliczenia były przeprowadzone dnie dla skalara x(1),x(2),... ale dla wektorów x(:,1),x(:,2),... Jeśli zastosujemy wbudowany generator normalnych liczb losowych '''normrnd''' to możemy także skorzystać z możliwości wygenerowania wielu liczb za jednym wywolaniem np. '''normrnd (0,1,M,1)'''. Postępując w ten sposób linijka:
| |
- | <source lang='matlab'>
| |
- | x(:,i)=x(:,i-1) + sqrt(2*D*h)* normrnd (0,1,M,1);
| |
- | </source>
| |
- | przedstawia wekrorowy zapis jednego kroku całkowania M stoschastychnych równań różniczkowych. Nie tylko upraszcza to zapis, ale także przyśpiesza obliczenia w przypadku stosowania interpretowanego języka jakim jest GNU Octave lub Matlab.
| |
- | [[Plik:wiener_Pxt.png|thumb|360px|Proces Wienera]]
| |
- | Cały program może wyglądać tak:
| |
- | <source lang='matlab'>
| |
- | clear all
| |
- | close all
| |
- | N=100;
| |
- | M=20000;
| |
- | T=4;
| |
- | h=T/N;
| |
- | clear x
| |
- | x=zeros(M,N);
| |
- | D=.81;
| |
- | for i=2:N
| |
- | x(:,i)=x(:,i-1) + sqrt(2*D*h)* normrnd (0,1,M,1);
| |
- | endfor
| |
- | </source>
| |
- |
| |
- | a procedura rysująca wykres histogramu położeń M cząstek po czasie t i porównująca
| |
- | ja z odpowiednim rozkładem Gaussa:
| |
- | <source lang='matlab'>
| |
- | n=10
| |
- | t=n/N*T
| |
- | xmax=4;
| |
- | h1=.2;
| |
- | hist(x(:,n),[-xmax:h1:xmax],1/h1)
| |
- | hold on
| |
- | fplot(@(xx) normpdf(xx,0,sqrt(2*D*t)),[-xmax,xmax],200,'ro-')
| |
- | hold off
| |
- | </source>
| |
- |
| |
- | Czyli w czasie funkcja rozkładu jest rozpływającym się po całej prostej rozkładem Gaussa. Poniższa procedura rysuje rozkłady w czterech momentach czasu:
| |
- | [[Plik:wiener_Pxt3.png|thumb|360px|Rozkład P(x,t) dla Procesu Wienera w czterech następujących po sobie chwilach. Widać rozwywanie się piku Gaussowskiego.]]
| |
- | <source lang='matlab'>
| |
- | for idx=1:4
| |
- | n=5+(idx-1)*20
| |
- | t=n/N*T
| |
- | subplot(2,2,idx)
| |
- | xmax=10;
| |
- | h1=.2;
| |
- | hist(x(:,n),[-xmax:h1:xmax],1/h1)
| |
- | hold on
| |
- | fplot(@(xx) normpdf(xx,0,sqrt(2*D*t)),[-24,24],200,'ro-')
| |
- | hold off
| |
- | legend ( sprintf("t=%2.0f",t) )
| |
- | endfor
| |
- | </source>
| |
- |
| |
- |
| |
- | === Niesymetryczny Proces Wienera: dyfuzja ze stałym dryftem ===
| |
- |
| |
- | Dyfuzja z dryfem jest granicznym przypadkiem [[PIZL:Proces_Wienera_i_proces_dyfuzji#Przypadek_niesymetryczny:_dyfuzja_z_dryfem|niesymetryczngo błądzenia przypadkowego]].
| |
- |
| |
- | Równanie dyfuzji z dryfem ma postać:
| |
- |
| |
- | <math>\frac{\partial P(x, t)}{\partial t} = -V \frac{\partial P(x, t)}{\partial x} + D \frac{\partial^2 P(x, t)}{\partial x^2} </math>
| |
- |
| |
- | i jest ono równoważne stochstycznemu równaniu różniczkowemu:
| |
- | <math>dX(t)= Vdt +\sqrt{2 D} dW(t)\;</math>.
| |
- |
| |
- | Łatwo możemy zmodyfikować program symulujący proces symmetryczny na niesymetryczny:
| |
- | [[Plik:wiener_dryf.png|thumb|360px|Niesymetryczny proces Wienera:dyfuzja z dryfem. Zielonym kolorem jest zaznaczona prosta x=Vt, która jest dryfem procesu]]
| |
- | <source lang='matlab'>
| |
- | close all
| |
- | clear all
| |
- | h=0.01;
| |
- | N=1200;
| |
- | x(1)=0;
| |
- | D=0.9;
| |
- | V=1.2;
| |
- | for i=2:N
| |
- | x(i)=x(i-1)+ V*h +sqrt(2*D*h)* normrnd (0,1,1,1);
| |
- | endfor
| |
- | plot((1:N)*h,x,'-',(1:N)*h,(1:N)*h*V,'-')
| |
- | </source>
| |
- |
| |
- | === Dyfuzja z dryftem: rozkład P(x,t) ===
| |
- |
| |
- | [[Plik:wiener_Pxt_dryf.png|thumb|360px|Rozkład P(x,t) dla Procesu Wienera w czterech następujących po sobie chwilach t = 0.2,1,1.8,2.6. Widać rozmywanie się piku Gaussowskiego.]]
| |
- | Podobnie jak w przypadku symetrycznego procesu Wienera rozważmy cząstkę starującą w punkcie <math>x=0</math> w czasie <math>t=0</math>. W języku gęstości prawdopodobieństwa oznacza to
| |
- |
| |
- | <math>P(x, 0) = \delta(x)\;</math>
| |
- |
| |
- | Jego rozwiązaniem na równania dyfuzji z dryfem prostej z takim warunkiem początkowym jest rozkład Gaussa:
| |
- |
| |
- | <math>P(x, t) = \frac{1}{\sqrt{4\pi D t} }\; \exp \left[ - \frac{(x-Vt)^2}{4Dt}\right],</math>
| |
- | który ma zależne od czasu odchylenie standardowe jak i wartość średnią.
| |
- |
| |
- | <source lang='matlab'>
| |
- | clear all
| |
- | close all
| |
- | N=100;
| |
- | M=20000;
| |
- | T=4;
| |
- | h=T/N;
| |
- | clear x
| |
- | x=zeros(M,N);
| |
- | D=.81;
| |
- | V=5.0
| |
- | for i=2:N
| |
- | x(:,i)=x(:,i-1) + V*h + sqrt(2*D*h)* normrnd (0,1,M,1);
| |
- | endfor
| |
- | hold on
| |
- | for idx=1:4
| |
- | n=5+(idx-1)*20
| |
- | t=n/N*T
| |
- | xmax=20;
| |
- | h1=.2;
| |
- | hist(x(:,n),[-2:h1:xmax],1/h1)
| |
- | fplot(@(xx) normpdf(xx,V*t,sqrt(2*D*t)),[-2,xmax],200,'ro-')
| |
- | endfor
| |
- | hold off
| |
- | </source>
| |
- |
| |
- | ===Proces Ornsteina Uhlenbecka===
| |
- |
| |
- | Jest to proces który spełnia równanie stochastyczne:
| |
- | [[Plik:OU20.png|thumb|360px|Dwadzieścia trajektorii procesu Ornsteina Uhlenbecka dla parametrów D=.01,k=1.0 i <math>\mu=1.0</math>. Dla czasów 0..2 widoczna jest relaksacja do stanu stacjonarnego, dla większych czasów fluktuacje wokół tego stanu.]]
| |
- |
| |
- | <math>dX(t)= -k ( X(t)-\mu ) dt + \sqrt{2 D } dW(t)\;</math>.
| |
- |
| |
- | Jest to proces z linowym dryfem. Posiada on stan stacjonarny (w przeciwieństwie do procesu Wienera). Jest modelem ciała zawieszonego sprężyście poddanego równowagowym fluktuacjom termicznym. Parametr <math>\mu</math> reprezentuje równowagową (stacjonarną) wartośc średnią procesu a k jest szybkością relaksacji do stanu równowagi.
| |
- |
| |
- | <source lang='matlab'>
| |
- | clear all
| |
- | close all
| |
- | N=500;
| |
- | M=20;
| |
- | T=14;
| |
- | h=T/N;
| |
- | clear x
| |
- | x=zeros(M,N);
| |
- | D=.01;
| |
- | k=1.0;
| |
- | mu=1.0;
| |
- | for i=2:N
| |
- | x(:,i)=x(:,i-1) - k*(x(:,i-1)-mu)* h + sqrt(2*D*h)* normrnd (0,1,M,1);
| |
- | endfor
| |
- | plot((1:N)*h,x,'r-')
| |
- | </source>
| |
Wersja z 12:42, 7 cze 2010
Stochastyczne równania różniczkowe
W tym rozdziale zostaną opisane metody numeryczne, które służa do rozwiązywania stochastycznych równań różniczkowych typu:
\(\frac{dX(t)}{dt} = F(X(t), t) + G(X(t), t)\Gamma(t)\)
gdzie F i G to dowolne funkcje, a \(\Gamma(t)\) jest procesem losowym. Najczęstszym przypadek to taki w którym \(\Gamma(t)\) to biały szum Gaussowski. Tak zapisane równanie nie jest precyzyjnie określone ze względu na dylemat Stratonowicza-Ito. Dlatego poprawne jest zapisanie równanie Ito w postaci:
\(dX(t)= F(X(t), t)dt + G(X(t), t) dW(t)\;\)
Nie zmienia to ogólności, gdyż jak wiadomo każde równanie zapisane w interpretacji Stratonowicza ma swój odpowiednik Ito. Dla potrzeb metod numerycznych będziemy rozpatrywać zawsze równania Ito, a jeśli pojawią się równania Stratonowicza to będziemy je transpormować do postaci Ito.
Schemat Eulera dla równań stochastycznych
Najprostszą metodą aproksymacji numerycznej równania stochastycznego jest podobnie jak w przypadku równań różniczkowych zwyczajnych jest schemat Eulera. Część deterministyczną równania stochastycznego traktujemy w taki sam sposób jak w schemacie Eulera dla równań różniczkowych zwyczajnych. Niech h oznacza krok całkowania i oś czasowa będzie zdyskretyzowana na przedzialy \(t_{i-1},t_{i},t_{i+1}\) oraz \(h=t_{i}-t_{i-1}\). Wtedy część deterministyczna równania stochastycznego przyjmuje postać:
\(X(t_i) = X(t_{i-1}) + \int_{t_{i-1}}^{t_{i}} F(X(t), t) dt \simeq X(t_{i-1}) + F(X(t_{i-1}, t_{i-1}) h \)
Aby całkować część stochastyczną potrzebujemy formuły na przyrost skończony procesu Wienera:
\( \int_{t_{i-1}}^{t_{i}} \Gamma(t) dt =\;\int_{t_{i-1}}^{t_{i}} dW(t) = W(t_{i})-W(t_{i-1})\)
Wiemy, że proces Wienera jest procesem o przyrostach niezależnych, które są gaussowską zmienna losową o zerowej wartości średniej i wariancji \(2(t_{i} − t_{i-1})=2 h\).
Tak więc widać, że w schemacie Eulera całkę typu \( \int_{t_{i-1}}^{t_{i}} \Gamma(t) dt \) należy zastąpić w każdym kroku całkowania gaussowską zmienną losową o wariancji proporcjonalnej do kroku całkowania h. Ponieważ z reguły dysponujemy gaussowskich generatorem liczb losowych o jednostkowej wariancji N(0,1), można zapisać:
\( \int_{t_{i-1}}^{t_{i}} \Gamma(t) dt = \sqrt{2 h} N(0,1) \)
Ten zapis pokazuje też ważną cechę przy obliczaniu aproksymacji rozwiązań stochastycznych - najniższy rząd w h jest nie O(h) ale \(O(h^{1/2})\).
Ponadto z takiego sformułowania widać też, że zmiany procesu Wienera w stosunku do przyrostów czasu są rozbieżne w granicy \(h\to 0\).
Korzystając w powyższych faktów, możemy zapisać pełny schemat Eulera dla równania stochastycznego (Ito):
\(X(t_i) = X(t_{i-1}) + F(X(t_{i-1}, t_{i-1}) h + \sqrt{h}G(X(t_{t-1}), t_{t-1}) N(0,1)\).
Gaussowskie zmienne losowe możemy otrzymać np. korzystając z algorytmu Box-a Mullera.
Schemat Eulera dla układu równań stochastycznych
Schemat Eulera można uogólnić na układy równań stochastycznych. Niech \(\mathbf X(t)\) będzie wektorem o składowych \[\mathbf X(t)=( X^1(t),X^2(t),...,X^n(t))\]. Układ równań stochastycznych (Ito) można zapisać w ogólnej postaci:
\[d X^i(t)= F^i(\mathbf X(t), t)dt + \sum_{j=1}^{n} G^{i,j}(\mathbf X(t), t) dW^j(t)\;\]\; j=1,2,...,n,
gdzie \(W^i(t),\;W^j(t)\) są niezależnymi procesami Wienera dla \(i\neq j\), \(F^i\) oznacza wektor drytfu a \(G^{i,j}\) jest macierzą \(n \times n\) funkcji.
Wtedy schemat Eulera ma postać:
\( X^j(t_i) = X^j(t_{i-1}) + F^j (\mathbf X(t_{i-1}), t_{i-1}) h + \sqrt{h} \sum_{k=1}^{n} G^{j,k}(\mathbf X(t), t) N^k(0,1)\; \)