Zajęcia IV: Ergodyczność

Rozwiązania i komentarze do laboratoryjnych zajęć z Modelowania Matematycznego

Zad. 1. (Zastosowanie Tw. Birkhoffa)

Oczywistym jest, że  obrót na okręgu zachowuje miarę Lebesgue' a.
Co więcej z zajęć I i IV wykładu możemy wnioskować, że  obrót o kąt niewspółmierny do π jest przekształceniem ergodycznym.

I. definujemy funkcję zwracającą orbitę długości n startującą z punktu (Subscript[x, 0])∈(0,2π) w funkcji obrotu o kąt ϕ.

orbitaobrotu[x0_,ϕ_,n_]:=NestList[N[Mod[#+ϕ,2 Pi]]&,x0,n-1]

Losujemy jeszcze przykładowe warości (Subscript[x, 0]) oraz ϕ:

x0=RandomReal[{0,2Pi}];
ϕ=RandomReal[{0,2Pi}];

II. Za pomocą funkcji Map[] będziemy stosować zadaną funkcję do każdej współrzędnej wygenerowanej orbity.:

Przykładowo :

Map[ Cos,orbitaobrotu[x0,ϕ,5]]

Zapisujemy funkcję wyliczającą sumy z twierdzenia Birkhoffa pomnożone przez λ ([0, 2 π]) = 2 π

aproksymacja[fun_,x0_,ϕ_,n_]:=2 Pi Total[Map[fun,orbitaobrotu[x0,ϕ,n]]]/n

II Obliczamy zadane całki

a)  (∫_0^(2π))sin(x)dx

Wartość dokładna

Integrate[Sin[x],{x,0,2Pi}]

wartości przybliżone

Table[a=aproksymacja[Sin,x0,ϕ,10^i];  "Dla n="<>ToString[Superscript[10,ToString[i]],StandardForm]<>" wartość przybliżona to: "<>ToString[a,TraditionalForm]<>" błąd bezwzględny wynosi:"<>ToString[N[Abs[a]],TraditionalForm],{i,0,4}]

a)  (∫_0^(2π))(cos^2)(x)dx

f2[x_]:=(Cos[x])^2

Wartość dokładna

N[Integrate[f2[x],{x,0,2Pi}]]

wartości przybliżone

Table[a=aproksymacja[f2,x0,ϕ,10^i];  "Dla n="<>ToString[Superscript[10,ToString[i]],StandardForm]<>" wartość przybliżona to: "<>ToString[a,TraditionalForm]<>" błąd bezwzględny wynosi:"<>ToString[N[Abs[a-Pi]],TraditionalForm],{i,0,4}]

Zad. 2. (rodzina kwadratowa raz jeszcze)

Nasza stara znajoma, funkcja z rodziny kwadratowej:

f=#1 #2(1-#2)&

Losujemy wektor m liczb pseudolosowych z [0, 1]. Tu dla przykładu m=4.

RandomReal[1,4]

Działamy na każdej współrzędnej tego wektora n=10 razy naszą funkcją

dane[λ_,n_,k_] :=Nest[f[λ,#]&,RandomReal[1,n],k]
dane[4.,4,10]

Tworzymy histogram :

Histogram[dane[3.3,500,1],10,"PDF"]

Zobaczmy jak zachowuje się w miarę iteracji

seed=RandomReal[1,1000];
Animate[Histogram[Nest[f[4.,#]&,seed,n],50,"PDF",PlotRange->{0,4}],{n,0,10,1}]

Stwórzmy więc listę, żeby porównać ze spodziewaną gęstością

HistogramList[Nest[f[4.,#]&,seed,30],50,"PDF"]

Rysowanie otrzymanej listy za pomocą ListPlot[]

hl= HistogramList[Nest[f[4.,#]&,seed,30],50,"PDF"];
estymmiara=ListPlot[Partition[Riffle[hl[[1]],hl[[2]]],2],Joined->True,Mesh->All]

Gęstość niezmiennicza dla rodziny logistycznej: (1/(π(Sqrt[x(1-x)])))

Show[
Plot[1/(Pi*Sqrt[x*(1-x)]),{x,0,1},PlotStyle->Red,AxesOrigin->{0,0}],
estymmiara]

Sprawdźmy jeszcze czy wydłużenie trajektorii i próbki może pomóc

Show[
Plot[1/(Pi*Sqrt[x*(1-x)]),{x,0,1},PlotStyle->Red,AxesOrigin->{0,0},Frame->True,GridLines->Automatic],
ListPlot[Partition[Riffle[#[[1]],#[[2]]],2],Joined->True,Mesh->All]]&@HistogramList[Nest[f[4.,#]&,RandomReal[1,1000000],5],100,"PDF"]

Show[
Plot[1/(Pi*Sqrt[x*(1-x)]),{x,0,1},PlotStyle->Red,AxesOrigin->{0,0},Frame->True,GridLines->Automatic],
ListPlot[Partition[Riffle[#[[1]],#[[2]]],2],Joined->True,Mesh->All]]&@HistogramList[Nest[f[4.,#]&,RandomReal[1,1000],5000],10,"PDF"]

Zad. 3.

Gęstość niezmiennicza dla rodziny logistycznej z parametrem λ=4:
 h(x)=(1/(π(Sqrt[x(1-x)])))

a) Dokładna wartość oczekiwana dana przez tą miarę :

Integrate[Sqrt[x]*(1/(π*Sqrt[x*(1-x)])),{x,0,1}]
N@%

b) Poprzez iteracje

Sprawdźmy, że miara absolutnie ciągła h (x) jest probabilistyczna

Integrate[(1/(π*Sqrt[x*(1-x)])),{x,0,1}]

Funkcja zwracająca orbitę o długości n dla rodziny logistycznej

orbitalogistyczna[x0_,n_]:=NestList[N[4*#*(1-#)]&,x0,n-1]

Losujemy punkt początkowy (dla  p.w. elementów (0, 1) metoda powinna działać)

ppocz=RandomReal[];

Wyliczamy wartość (1/n)(∑_(i=0)^(n-1))f((S^i)(x)) dla dużych n, gdzie f(x)=(Sqrt[x]).

aproksymacjalog[x0_,n_]:=Total[Map[Sqrt,orbitalogistyczna[x0,n]]]/n

Table[q=aproksymacjalog[ppocz,10^i]; "Dla n="<>ToString[Superscript[10,ToString[i]],StandardForm]<> " wartość przybliżona to: "<>ToString[q]<>" błąd bezwzględny wynosi:"<>ToString[N[Abs[q-2/Pi]]],{i,0,7}]

Zad. 4.

Zauważmy najpierw, że takie układy dynamiczne można zapisac w postaci 0-1 macierzy stochastycznych. A zatem pytanie sprowadza się do wygenerowania wszystkich macierzy stochastycznych 3 3 . Zapytajmy najpierw: ile ich jest?

Każdy ze stanów może prowadzić do każdego z trzech stanów, a zatem

(3^3)

Zapiszmy to najpierw korzystając z notacji DirectedEdge  -> (ESC-de-ESC) , szukając wszelkich możliwych kombinacji

li=Flatten[Table[{1\[DirectedEdge]i,2\[DirectedEdge]j,3\[DirectedEdge]k},{i,3},{j,3},{k,3}],2](*Spłaszczamy listę poza najniższym poziomem*)
Length@li

Zauważmy, że funkcja Graph pomoże nam w wygodny sposób narysować grafy przejścia

Graph[#,VertexLabels->"Name"]&/@li

Proszę zwrócić uwagę, że elementy powyższej listy nie są obiektami graficznymi tylko grafami (zachęcam do sprawdzenia możliwości jakie daje prawoklik myszy).

Narysujmy to nieco ładniej

Graph[#,VertexShapeFunction->"Hexagon",VertexSize->0.2,VertexLabels->Placed["Name",Center],VertexLabelStyle->Directive[12,White,Bold,FontFamily->"Courier"],EdgeStyle->Directive[Red,Thick]]&/@li

Jak jednak otrzymać macierze sąsiedztwa tych grafów?

Najpierw uprośćmy maksymalnie naszą składnię:

Graph/@li

Teraz sprawdźmy działanie funkcji AdjacencyMatrix

AdjacencyMatrix/@(Graph/@li)

SparseArray to funkcja, która przechowuje w pamięci rzadkie macierze, zmieńmy je na bardziej czytelne:

MatrixForm/@(AdjacencyMatrix/@(Graph/@li))

Wygląda to bardzo ładnie, ale nie jest zbyt użyteczne, zamieńmy to więc na dobrze znane nam listy

macierze=Normal/@(AdjacencyMatrix/@(Graph/@li))

Sprawdźmy teraz do czego prowadzi iteracja każdej macierzy

MatrixPower[#,2]&/@macierze

Dyskwalifikujemy te macierze, których grafy nie są spójne w słabym sensie (Dlaczego?) w tym celu używamy funkcji AdjacencyGraph, która działa odwrotnie do AdjacencyMatrix, zwracając graf dla danej macierzy

Select[macierze,WeaklyConnectedGraphQ[AdjacencyGraph[#]]&]
AdjacencyGraph/@%
Length@%

Proszę zwrócić uwagę, że nie sprecyzowaliśmy jeszcze jaką miarę rozważamy na podanych zbiorach. Okazuje się, że jeśli oczekujemy, żeby miara każdego z wierzchołków była dodatnia to musimy wykluczyć wszystkie funkcje poza dwoma typu

A zatem mamy nasz ergodyczny w sensie dynamicznym układ (proszę sprawdzić, że z miarą {1/3,1/3,1/3} spełnia założenia) Czy jest on ergodyczny w sensie procesu Markowa?

Jest nieprzywiedlny!

Nie jest nieokresowy :(

Zad. 5.

Rozważmy błądzenie losowe po grafie pełnym

CompleteGraph[5]

am=Normal@AdjacencyMatrix[CompleteGraph[5]]

Potrzebujemy jeszcze unormować tę macierz, żeby była stochastyczna

m=#/Total[#]&/@am

Iterujemy macierz

Table[MatrixPower[m,n],{n,10}]//N

Z powyższego od razu widzimy, że rozkład ergodyczny to {0.2, 0.2, 0.2, 0.2, 0.2}, ale upewnijmy się jeszcze

Eigensystem[m]

Rzeczywiście, z wartością własną 1 wiąże się wektor własny {1, 1, 1, 1, 1}, który po unormowaniu równy będzie {1/5, 1/5, 1/5, 1/5, 1/5}