Metody komputerowe modelowania 
układów fizycznych

Wykład dla Studium Doktoranckiego Wydziału Matematyczno-Przyrodniczego 
Uniwersytetu Jana Kochanowskiego w Kielcach

2015/16

Wojciech Broniowski

Najnowszą wersję tego pliku znajdziesz na mojej stronie domowej http://www.ujk.edu.pl/~broniows

Wersja 16.1, 24 II 2016, Mathematica 10.3.1

(Dla raczkujących: poniżej kliknij dwa razy na klamry z haczykiem po prawej stronie aby otworzyć/zamknąć komórkę. Podobnie jest z trójkącikami po lewej stronie.)

Ustawienie ścieżek do czytania danych - poniżej zmień na swoją ścieżkę, by wczytywać dane.

(* Windows, przykładowo: *)
(* datapath="d:\\moj_katalog\\"; *)

(* unix *)
datapath="/home/broniows/local/Computer_Methods_2015/"

$Path=AppendTo[$Path,datapath];

Globalne opcje :

SetOptions[{Plot,LogPlot,ListPlot},{Frame->True,PlotStyle->Thick,FrameTicksStyle->Directive[FontSize->14,FontWeight->Bold],LabelStyle->Directive[FontSize->18,FontWeight->Bold]}];
SetOptions[{Plot3D,RegionPlot3D},{TicksStyle->Directive[FontSize->14,FontWeight->Bold],LabelStyle->Directive[FontSize->18,FontWeight->Bold]}];

Wykład 1 
Liczby pseudolosowe i całkowanie Monte Carlo

Wbudowany generator liczb pseudolosowych

Własności dobrego generatora liczb pseudolosowych:

* brak korelacji
* długi okres
* niezależność właściwości od wyboru ziarna
* szybkość

MLCG - Multiplicative Linear Congruential Generators

(Subscript[x, 0])=ziarno   - np. zależy od daty i czasu
(Subscript[x, n+1])=((Subscript[x, n])*a+b) mod c  - a, b, c odpowiednio dobrane liczby naturalne

r=(Subscript[x, n])/c daje liczbę z przedziału [0,1],   c bardzo duże, rzędu okresu

Przykładowo,

a=1103515245, b=12345, c=(2^31) (funkcja rand, UNIX 32 bit)
a=(13^13), b=0, c=(2^59) (Numerical Algorithms Group)

r:=Random[]

Dzieki dwukropkowi za każdym wywołaniem inna liczba!

r

Timing[t=Table[r,{10000000}];]

Rozkład jest istotnie jednorodny (nieodzowne fluktuacje związane są ze skończoną krotnością próbki.

Histogram[t,10,"PDF"]

Całkowanie metodą Monte Carlo

Funkcja f(x), dodatnio określona, zadana na przedziale [a,b], górne ograniczenie M

f[x_]=Sin[x^2]^4
M=1
a=0
b=Sqrt[Pi]
Plot[f[x],{x,a,b}]

Spróbujmy ją scałkować analitycznie

Integrate[f[x],{x,a,b}]
%//N

Wynik jest analityczny i zawiera funkcję Fresnela. Scałkujmy numerycznie:

NIntegrate[f[x],{x,a,b}]

Teraz obliczymy to samo “na piechotę”, tj. z pomocą metody Monte Carlo. W tym celu losujemy punkt z prostokąta [a,b]x[0,M]

po:={Random[Real,{a,b}],Random[Real,{0,M}]}

po

Jeśli punkt jest pod wykresem funkcji, to go akceptujemy (nin++), w przeciwnym razie odrzucamy. Wykonujemy np losowań. Wartość całki jest przybliżana liczbą nin/np  M (b-a).

np=1000000;
Module[{},nin=0;Do[{x,y}=po;If[y<f[x],nin++],{i,1,np}];nin/np  M (b-a)//N]

Porównując z wcześniejszym wynikiem, jest zupełnie nieźle. Jednak w metodzie Monte Carlo względny błąd skaluje sie jak 1/n, gdzie n jest liczbą losowań, zatem uzyskanie większej dokładności jest kosztowne. 

Poniżej wariant z zapamiętywaniem punktów celem wizualizacji

np=10000

Module[{x,y},nin=0;tt={};Do[{x,y}=po;AppendTo[tt,po];If[y<f[x],nin=nin+1],{i,1,np}];nin/np  M (b-a)//N]

plo=ListPlot[tt];

plof=Plot[f[x],{x,a,b},PlotStyle->{Thickness[.01],Red}];

Show[plo,plof,DisplayFunction->$DisplayFunction]

Igła Buffona (1733, rozw. 1777) - przykład Monte Carlo oraz grafiki

Quit[]

Słynny przykład “analogowego” rzucania igły o długości l na paper w linie odległe od siebie o d. Zliczamy stosunek przypadków, gdzie igła przetnie linię, do liczby wszystkich rzutów. Pozwala to znaleźć liczbę π!

L=1 (* długość igły *)
d=1 (* odległość między liniami *)

Plot3D[x+L/2 Cos[ϕ],{x,-d/2,d/2},{ϕ,-Pi,Pi},AxesLabel->{"x","ϕ"}]

Plot3D[{-d/2,x+L/2 Cos[ϕ],d/2},{x,-d/2,d/2},{ϕ,-Pi,Pi},AxesLabel->{"x","ϕ"},PlotStyle->Opacity[.7],Mesh->True]

Z symetrii, wystarczy rozważać górny i prawy obszar

Plot3D[{x+L/2 Cos[ϕ],d/2},{x,-d/2,d/2},{ϕ,0,Pi},AxesLabel->{"x","ϕ"},PlotStyle->Opacity[.6],Mesh->False]

Dalej bierzemy dowolne L, d może pozostać = 1 jako skala

Clear[L,x]

Solve[x+L/2 c == d/2,x]  (* warunek na linię rozgraniczającą, c=Cos[phi] *)

D[Cos[phi],phi] (* jakobian *)
jac=1/Sqrt[1-c^2]

i1=Integrate[1,{x,1/2(1-c L),1/2}]

Liczba rzutów z przecięciami, czynnik 2 z przecięć z dolną linią (patrz wykresy powyżej):

i2=2 Integrate[jac i1,{c,0,1},Assumptions->{l>0}]

i10=Integrate[1,{x,-1/2,1/2},{x,0,1}]

Liczba wszystkich rzutów :

i20=Integrate[jac i10,{c,0,1}]

Z tego stosunku dostajemy π:

i20/i2  2L/d

tj. dla dowolnych d i L otzymujemy π = 2 L/d  (liczba wszystkich rzutów)/ (liczba rzutów z przecięciami)

Teraz to samo metodą Monte Carlo

phi:=Random[Real,{-Pi/2,Pi/2}]
x:=Random[Real,{-10d,10d}]
y:=Random[Real,{-10d,10d}]

Igła o kącie phi i położeniu środka (x,y)

igla:={phi,x,y}

igla

Warunek przecięcia linii (prosta geometria, z której to wynika, jest pozostawiona jako ćwiczenie.

tnie[ig_]:=Module[{m},m=Mod[ig[[3]],d];m=Min[m,1-m];l/2 Cos[ig[[1]]]>m]

tnie[igla]

gr[ig_]:={If[tnie[ig],Red,Blue],Line[{{ig[[2]]+l/2 Sin[ig[[1]]],ig[[3]]+l/2 Cos[ig[[1]]]},{ig[[2]]-l/2 Sin[ig[[1]]],ig[[3]]-l/2 Cos[ig[[1]]]}}]}

Show[Graphics[gr[igla]],ImageSize->100]

igly={}

Do[igly=AppendTo[igly,gr[igla]],{i,1,1000}]

tt=Table[-10+i,{i,0,20}]

canvas=Plot[tt d,{x,-10d,10d},AspectRatio->1,Frame->False,Axes->False]

Show[canvas,Graphics[igly],Axes->False]

igly={}

ListAnimate[Table[Show[canvas,Graphics[igly=AppendTo[igly,gr[igla]]],Axes->False],{i,1,100}]]

... i wreszcie obliczenie π w n rzutach

s=0 (* liczba rzytów z przecięciami *)
n=4000000 (* liczba wszystkich rzutów *)
Timing[Do[If[tnie[igla],s++],{i,1,n}];]
res = 2 l n/(d s)//N
100 (res-Pi)/Pi "%"

Wbudowane rozkłady

Quit[]

RandomVariate[NormalDistribution[0,4]]

Show[Histogram[Table[RandomVariate[NormalDistribution[0,1]],{20000}],50,"PDF"],Plot[Exp[-x^2/2]/Sqrt[2Pi],{x,-4,4}]]

Dowolny rozkład - odwrotna dystrybuanta

Zobaczmy, jak transformują się rozkłady prawdopodobieństwa przy zamianie zmiennej. Mamy
p(y) dy  =f(x)dx    
(ale f(x)=1, bo x ma z założenia rozkład jednorodnym). Zatem p(y) dy  =dx oraz dystrybuanta zmiennej y wynosi
P(y)=x
oraz
y=(P^(-1))(x)

A zatem aby otrzymać dowolny rozkład należy losować zmienną jednorodną x, a następnie "obłożyć ją" odwrotną dystrybuantą żądanego rozkładu.

Quit[]

Przykład :

f[x_]=x^2

P[y_]=Integrate[f[x],{x,0,y},Assumptions->y>0]

Pinv[z_]=y/.Solve[P[y]==z,y][[2]]

gd:=Random[]^(1/3)

Show[Histogram[Table[gd,{100000}],1000,"PDF"],Plot[3x^2,{x,0,1},PlotStyle->Red],PlotRange->All]

Quit[]

?f

f[x_]=1/(1+x^2)
a=-1
b=1
M=1

x

no=Integrate[f[xx],{xx,a,b}]

pp=Plot[f[x]/no,{x,-1,1},PlotRange->{0,1}]

ra:=Module[{},y=2;x=0;While[y>f[x],x=Random[Real,{a,b}];y=Random[Real,{0,1}]];x]

ra

T=Table[ra,{1000}];
Show[Histogram[T,10,"PDF"],pp]

Dowolny rozkład - von Neumann

Quit[]

Rozkład sin(x)

f[x_]=Sin[x]/2;
low=0;
up=Pi;
maxim=1;

ns:=Module[{x,y},y=∞;x=low;While[y>f[x],{x,y}={RandomReal[{low,up}],RandomReal[{0,maxim}]}];x]

ta2=Table[ns,{i,1,10000}];

plh=Histogram[ta2,50,"PDF"];
plf=Plot[f[x],{x,low,up},PlotStyle->{Thickness[0.005]}];
Show[plh,plf]

Chi kwadrat i fitowanie danych do modelu

n=15;

x:=Sum[RandomVariate[NormalDistribution[0,1]]^2,{n}]
y:=RandomVariate[ChiSquareDistribution[n]]

tx=Table[x,{100000}];
ty=Table[y,{100000}];
Histogram[{tx,ty},{1},"PDF"]

Quit[]

<<"ErrorBarPlots`";

myst=OpenRead["PH200s.txt"]

(* ustaw na początku *)
SetStreamPosition[myst,0]

(* opuść 6 linii *)
Skip[myst,Record,6]

tt={}

Do[tt=AppendTo[tt,Read[myst,{Real,Word,Real,Word,Real,Word,Real,Word,Real,Word,Real,Word,Real,Word,Real,Word,Real}]],{i,1,28}]

Tak wygląda wczytana tabela

tt //MatrixForm

"Wyławianie" istotnych kolumn

t1=Table[{tt[[i]][[1]],tt[[i]][[3]],tt[[i]][[5]]},{i,1,Length[tt]}];

t1//MatrixForm

plexp=ErrorListPlot[t1,PlotRange->All]

To samo dla wykresu logarytmicznego

t1plog=Table[{t1[[i]][[1]],Log[10,t1[[i]][[2]]],t1[[i]][[3]]/Abs[t1[[i]][[2]]]},{i,1,Length[t1]}];

plexppl=ErrorListPlot[t1plog,PlotRange->All,PlotStyle->PointSize[0]]

t2=Table[{tt[[i]][[1]],Log[10,tt[[i]][[3]]]},{i,1,Length[tt]}]

Clear[x]

Hipoteza modelowa

f[x_,A_,a_,b_,c_]=A (1+a x+b x^2 +c x^3)

opt=FindFit[t2,f[x,A,a,b,c],{A,a,b,c},x]

ff[x_]=f[x,A,a,b,c]/.opt

Show[Plot[ff[x],{x,0,3}],ListPlot[t2]]

(* χ kwadrat na stopień swobody *)

f[x_,A_,a_,b_]=A -a x -b x^2

chi2[A_,a_,b_]=Sum[(f[t1plog[[i]][[1]],A,a,b]-t1plog[[i]][[2]])^2/(t1plog[[i]][[3]])^2,{i,1,Length[t1]}]/(Length[t1]-3);

Minimalizacja chi kwadrat

opti1=FindMinimum[chi2[A,a,b],{A,10,20},{a,1.5,1.7},{b,.2,.3}]

aopt=a/.opti1[[2]]
bopt=b/.opti1[[2]]
Aopt=A/.opti1[[2]]
chimin=opti1[[1]]

pl=Show[Plot[f[x,Aopt,aopt,bopt],{x,0,3},PlotStyle->RGBColor[1,0,0]],ErrorListPlot[t1plog]]

ContourPlot[chi2[Aopt,a,b],{a,.9aopt,1.1 aopt},{b,.5bopt,1.5 bopt},ContourShading->False,PlotPoints->100,AxesLabel->{"a","b"},Contours->{chimin+3.53}]

RegionPlot3D[chi2[A,a,b]< chimin+3.53,{A,.8Aopt,1.2Aopt},{a,.8aopt,1.2 aopt},{b,.7bopt,1.3 bopt},AxesLabel->{"A","a","b"},PlotRange->Automatic, PlotPoints->60, Mesh->None,PlotStyle->{Red,Opacity[.3]}]

f[x,Aopt,aopt,bopt]

wei=Table[1/t1plog[[i,3]]^2,{i,1,Length[t1plog]}];

myfit=NonlinearModelFit[t2,A-a x - b x^2,{A,a,b},x,Weights->wei]

myfit["CorrelationMatrix"]//MatrixForm

myfit["ANOVATable"];
myfit["ParameterTable"]
pp=myfit["MeanPredictionBands"];
Show[Plot[pp,{x,0,3},Filling->{1->{2}},PlotStyle->{Thickness[0.00001]}],plexppl]
myfit["CorrelationMatrix"]//MatrixForm

ZADANIE 1

0) [Opcjonalne] Zainstalować Mathematikę i zrobić piersze rozdziały z zielonej księgi -- przejść mój przewodnik anty_7.nb (zob. moją stronę www) lub dowolny inny w miarę zaawansowany tutorial.
1) [obowiązkowe na za tydzień] Zaprogramuj wybrany multiplikatywny generator liczb pseudolosowych. 

************************************************

2) [onzt] Skonstruuj metodą a) odwrotnej dystrybuanty oraz b) von Neumanna generator rozkładu radialnego Woodsa-Saxona, f[r]=r^2/(1+Exp[r]). Porównaj, który jest szybszy. Wskazówka: w części a) należy stablicować odwrotną dystrybuantę i skonstruować jej interpolację z pomocą “Interpolation”.

Rozwiązanie 1

Quit[]

mySeed[s_]:=Module[{},xran=s;aran=1103515245;bran=12345;cran=2^31;Print["myRan seeded with ",s]]

myRan:=Module[{},xran=Mod[xran aran+bran,cran];xran/cran //N]

mySeed[33]

myRan

t=Table[myRan,{1000000}];
Histogram[t,10,"PDF"]

Rozwiązanie 2 a)

Quit[]

Teatr :

f[r_]=4Pi r^2/(1+Exp[(r-R0)/a])

norm=Integrate[f[r],{r,0,Infinity},Assumptions->a>0]

fn[r_]=f[r]/norm

pl=Plot[fn[r]/.{R0->1,a->.1},{r,0,2}]

F[y_,R0_,a_]=Integrate[fn[r],{r,0,y},GenerateConditions->False]

Plot[{F[y,1,.1],1},{y,0,2}]

Finv[x_]:=y/.FindRoot[F[y,1,.1]==x,{y,.5}]

Plot[Finv[x],{x,.1,.99}]

tt=Table[{x,Finv[x]},{x,0,.9999,.9999/1000}]

Finvint=Interpolation[tt]

Plot[{Finvint[x]/Finv[x]},{x,0,.999},PlotRange->{.9,1.1}]

w:=Finvint[Random[]]
Timing[tab=Table[w,{10000}];]

Show[Histogram[tab,30,"PDF"],pl]

Wykład 2  
Model Isinga i algorytm Metropolisa

Poradnik

Minikompendium fizyki statystycznej:

Układ jest z założenia w równowadze termicznej z rezerwuarem ciepła, h - zewnętrzne pole magnetyczne, M - magnetyzacja

Prawdopodobieństwo obsadzenia stanu i o energii (Subscript[E, i]) oraz magnetyzacji (Subscript[M, i]): (Subscript[P, i])=1/Z (e^(-((Subscript[E, i])-(Subscript[M, i]) h)/((Subscript[k, B])T)))
Normalizacja (funkcja rozdziału): Z= (∑_i)(e^(-((Subscript[E, i])-(Subscript[M, i]) h)/((Subscript[k, B])T))) - zawiera całą informację o układzie
β=1/T, stała Boltzmanna (Subscript[k, B])=1

średnia energia:<E>=(∑_i) (Subscript[P, i]) (Subscript[E, i])= - (∂/(∂β ))log Z

ciepło właściwe: (Subscript[C, V]) ≡ ((∂<E>)/(∂T)) = (β^2) [(1/Z) (((∂^2)Z)/(∂(β^2))) -  (1/(Z^2)) (((((∂Z)/(∂β))))^2) ]= (β^2) [<(E^2)> - <E((>)^2)] - tw. o fluktuacji i dysypacji (!)

średnia magnetyzacja: <M> = (∑_i) (Subscript[P, i]) (Subscript[M, i])=  (∂/(β∂h ))log Z

podatność magnetyczna: χ ≡ ((∂<M>)/(∂h))=β [<(M^2)> - <M((>)^2)] 


Magnetyki jako przykład 

* Prosty model magnetyzmu. 
* Przykład symulacji MC w fizyce statystycznej.
* Zastosowania w przejściach fazowych, zjawiskach krytycznych, kwantowej teorii pola (siatki). 

Magnetyki podlegają przejściom fazowym. Zachodzą rozbieżności pewnych wielkości (podatność magnetyczna) w temperaturze Curiego.
	niska temperatura ⟶ porządek ⟶ wspólny kierunek spinów ⟶ wysoka magnetyzacja
	wysoka temperatura ⟶ nieporządek ⟶ losowy kierunek spinów ⟶ brak magnetyzacji

1D model Isinga (Wihelm Lenz 1920, Ernst Ising 1925 )

Energia układu: E=-J  (∑_(<i,j>))(Subscript[S, i])(Subscript[S, j])
(Subscript[S, i])= +1 lub -1   - spin w węźle i siatki
<i,j>  -  najbliźsi sąsiedzi
J   - stała oddziaływania 

Jeśli J > 0  - ferromagnetyk, jeśli J < 0 - antiferromagnetyk (rozważymy przypadek ferromagnetyka).

Ising: brak przejścia fazowego w 1D.
Lars Onsager (1944): dokładne rozwiązanie 2D, występowanie przejścia fazowego jak w rzeczywistych ferromagnetykach.
3D - nierozwiązywalny analitycznie, symulacje.

Fizyka statystyczna

Magnetyzacja jest średnią termodynamiczną po stanach C, 

<M(T)>=(∑_C)p(C,T)M(C)

T - temperatura
C - konfiguracja
p(C,T) - prawdopodobieństwo konfiguracji C
M(C) - magnetyzacja konfiguracji C

Boltzmann ((Subscript[k, B])=1):

p(C)=(Z^(-1)) (e^(-E(C)/T))
Z=(∑_C) (e^(-E(C)/T)) - funkcja rozdziału
E(C) - energia konfiguracji C

Metoda Monte Carlo 

Naiwna metoda oparta o niezależne próbkowanie

<M>=(Underscript[(∑^n), i=1]) M((Subscript[C, i]))(e^(-E((Subscript[C, i]))/T))/(Underscript[(∑^n), i=1]) (e^(-E((Subscript[C, i]))/T))

nie zadziała! Większość losowo wybranych konfiguracji będzie mieć znikomą wagę, tj. nie natrafimy na chybił trafił na konfiguracje “dobre”

Metoda brutalna

(Subscript[N, A ])=6(x10^23) - tyle spinów mamy w układzie makroskopowym. Typowa siatka 2D: 24x24.

Liczba konfiguracji = (2^(24x24)) ≈ (10^172) 

(10^(9 ))konfiguracji/procesor/s x(10^(10 ))procesorów x (10^14)s/rok x (10^10)lat = (10^43) - tylko tyle sprawdzimy!

Próbkowanie istotności (importance sampling)

Gdybyśmy próbkowali konfiguracje o rozkładzie Boltzmanna, wtedy średnia Monte Carlo wyniosłaby

<M>=(Underscript[(∑^n), i=1]) M((Subscript[C, i]))/n

Nietety, rozkładu tego (czy też funkcji rozkładu) nie mamy - dopiero próbujemy go policzyć. 

Procesy Markowa

P(A,t) - prawdopodobieństwo konfiguracji A po “czasie” (liczbie kroków) t
W(A→B) - prawdopodobieństwo przejścia na jednostkę czasu ze stanu A do stanu B 

Równanie zysku/straty:

P(A,t+1)=P(A,t)+(∑_B)[W(B→A)P(B,t)-W(A→B)P(A,t)]  (**)

Dla dużych t zachodzi P(A,t)→p(A) - utrata pamięci warunków początkowych

Równowaga szczegółowa
 
Warunkiem wystarczającym dla (**) jest zatem

W(B→A)p(B)=W(A→B)p(A)

tj.

W(B→A)/W(A→B)=p(A)/p(B)=(e^(-ΔE/T))  (*)

ΔE=E(A)-E(B)

Algorytm Metropolisa (Metropolis, Teller, Rosenbluth, 1953)

Jak zadośćuczynić (*)? Jedna (z nieskończenie wielu) możliwości to

W(A→B)=(e^(-(E(B)-E(A))/T))        gdy  E(B)>E(A)  (wzbudza się z prawdopodobieństwem malejącym z energią wzbudzenia)
                   =1                                    gdy  E(B)⩽E(A) (przechodzi do stanu o niższej energii z prawdopodobieństwem 1) 

Wówczas jeśli E(B)>E(A), to

W(B→A)/W(A→B)=1/(e^((E(A)-E(B))/T))=(e^(-(E(A)-E(B))/T))= (e^(-ΔE/T)),

natomiast w przypadku E(B)<E(A) 

W(B→A)/W(A→B)=(e^((E(B)-E(A))/T))/1=(e^(-ΔE/T)),

czyli również w porządku!

Ogólne wymogi

Algorytm Monte Carlo musi spełniać

* równowagę szczegółową - spełniona dla algorytmu Metropolisa
* ergodyczmność (każda konfiguracja jest możliwa do uzyskania z dowolnej innej w skończonej liczbie kroków) - spełnienie zależy od modelu, temperatury, ogólnie jest to nietrywialny problem

Teatr

Quit[]

n=100 (* number od lattice sites *)
J=1 (* J=1 - ferro *)

init:=Module[{},s=Table[2RandomInteger[1]-1,{n}];]

dis:=ListPlot[s]

ras:=Show[Graphics[Raster[{(s+1.1)/3}]],ImageSize->800]

Energy and magnetization

ener:=-J(Sum[s[[i]]s[[i+1]],{i,1,n-1}]+s[[n]]s[[1]]) //N
mag:=Sum[s[[i]],{i,1,n}] //N

init
ras
ener
mag

Metropolis

flip[i_]:=s[[i]]=-s[[i]];

de[i_]:=2J s[[i]](s[[i-1]]+s[[i+1]])/; 1<i<n 
de[1]:=2J s[[1]](s[[n]]+s[[2]])
de[n]:=2J s[[n]](s[[n-1]]+s[[1]])

Test

(* test *)

ir=RandomInteger[{1,n}]

de[ir]
e1=ener;
flip[ir]
e2=ener;
e2-e1

(* Metropolis *)
init;
met[i_]:=Module[{},If[de[i]<0,flip[i]];]
sweep:=Module[{j},j=RandomInteger[{2,n}];Do[met[i],{i,j,n}];Do[met[i],{i,1,j-1}];]

met[i_,T_]:=Module[{},If[de[i]<0,flip[i],If[Exp[-de[i]/T]>Random[],flip[i]]]]
sw[T_]:=Module[{j},j=RandomInteger[{2,n}];Do[met[i,T],{i,j,n}];Do[met[i,T],{i,1,j-1}];]

T0=.2
nsweeps=1000
init;
taben={ener}
tabs={mag}
Do[sw[T0];AppendTo[taben,ener];;AppendTo[tabs,mag],{nsweeps}]
ListPlot[{taben,tabs},PlotRange->All]

ninit=1000; (* czas dochodzenia do rownowagi *)
ncoll=5000; (* przemiecenia pomiarowe *)

CV[T_]:=Module[{},Do[sw[T],{k,1,ninit}];scoll=0;scoll2=0;Do[sw[T];scoll=scoll+ener;scoll2=scoll2+ener^2,{k,1,ncoll}];(scoll2/ncoll-scoll^2/ncoll^2)/T^2]

chi[T_]:=Module[{},Do[sw[T],{k,1,ninit}];scoll=0;scoll2=0;Do[sw[T];scoll=scoll+mag;scoll2=scoll2+mag^2,{k,1,ncoll}];(scoll2/ncoll-scoll^2/ncoll^2)/T]

Table[CV[T],{T,.1,1,.1}]

ListPlot[%]

Z domu

Quit[]

(* liczba węzłów sieci *)
n=100

(* siła oddziaływania, +1 ferro, -1 ferri, można przyjąć jako skalę, więc |J|=1 *)
J=-1

(* losowa inicjalizacja tablicy *)
init:=t=Table[1-2Random[Integer],{i,1,n}]

init

(* sprawdzenie periodycznego "zawijania" - wygodne, aby poniżej pisać sumy w sposób jednolity *)
Table[{i,Mod[i+1,n,1]},{i,1,n}]

(* spin, energia i korelacja spinowa na węzeł *)
spin:=Sum[t[[i]],{i,1,n}]/n //N
energy:=-J Sum[t[[i]]t[[Mod[i+1, n,1]]],{i,1,n}]/n //N
gamma[m_]:=Sum[t[[i]]t[[Mod[i+m, n,1]]],{i,1,n}]/n (* - spin^2 *)//N

init;
spin
energy
(* {gamma[1],gamma[2],gamma[3],gamma[4]} *)

(* operacja przerzucenia spinu w węźle i *)
flip[i_]:=t[[i]]=-t[[i]]

(* zmiana energii przy przerzuceniu spinu w węźle i *)
del[i_]:=2J t[[i]](t[[Mod[i-1, n,1]]]+t[[Mod[i+1, n,1]]])

(* przeskok Metropolisa w węźle i przy temperaturze T *)
step[i_,T_]:=Module[{},If[del[i]<0,flip[i],If[Exp[-del[i]/T]>Random[],flip[i]]]]
(* jedno "przemiecenie", począwszy od losowo wybranego węzła r *)
sweepT[T_]:=Module[{r},r=Random[Integer,{1,n}];Do[step[i,T],{i,r,n}];Do[step[i,T],{i,1,r-1}];]

(* losowe wykonanie n przerzuceń - nie używane *)
(* sweepTran[T_]:=Do[step[Random[Integer,{1,n}],T],{n}] *)

Zobaczmy, jak to działa

init

energy

sweepT[.1]
energy

(* osiąganie równowagi termodynamicznej przy temperaturze T, liczba przemieceń = steps*)
etrace[T_,steps_]:=Module[{},ett={};stt={};Do[sweepT[T];AppendTo[ett,energy];AppendTo[stt,spin] (*; Print[energy, " ",spin],*),{k,1,steps}];Show[ListPlot[{Legended[ett,"energia/n"],Legended[stt,"spin/n"]}, PlotLabel->StringJoin["T=",ToString[T]],AxesLabel->{StyleForm["krok",FontSize->18],StyleForm["E/n, M/n",FontSize->18]},ImageSize->600 ,Joined->True,Frame->True,PlotRange->{-1.1,1.1}](*,Plot[Table[-1+4i/100,{i,0,50}],{x,0,steps}]*)]]

(* dodatkowo, możemy wizualizować losowy stan początkowy i stan po wielu przemieceniach *)
init;
Show[Graphics[Raster[{(t+1.2)/3}]],ImageSize->800,Background->LightBlue]
etrace[.95, 1000]
Show[Graphics[Raster[{(t+1.2)/3}]],ImageSize->800,Background->LightBlue]

Symulacje

ninit=500; (* liczba przemieceń do osiągnięcia równowagi *)
ncoll=2000; (* i nastepująca po niej liczba przemieceń pomiarowych *)

(* średni (po ncoll stanach w równowadze) spin, energia, ciepło właściwe i podatność magnetyczna *)

sav[T_]:=Module[{},Do[sweepT[T],{k,1,ninit}];scoll=0;Do[sweepT[T];scoll=scoll+spin,{k,1,ncoll}];scoll/ncoll]

eav[T_]:=Module[{},Do[sweepT[T],{k,1,ninit}];scoll=0;Do[sweepT[T];scoll=scoll+energy,{k,1,ncoll}];scoll/ncoll]

CV[T_]:=Module[{},Do[sweepT[T],{k,1,ninit}];scoll=0;scoll2=0;Do[sweepT[T];scoll=scoll+energy;scoll2=scoll2+energy^2,{k,1,ncoll}];(scoll2/ncoll-scoll^2/ncoll^2)/T^2]

chi[T_]:=Module[{},Do[sweepT[T],{k,1,ninit}];scoll=0;scoll2=0;Do[sweepT[T];scoll=scoll+spin;scoll2=scoll2+spin^2,{k,1,ncoll}];(scoll2/ncoll-scoll^2/ncoll^2)/T]

(* tablicowanie i wykres *)
teav=Table[{T,eav[T]},{T,.2,5,.2}];
tsav=Table[{T,sav[T]},{T,.05,5,.2}];

plmc=Show[ListPlot[teav,PlotStyle->Red,Joined->True,Frame->True],ListPlot[tsav,Joined->True],ImageSize->600,FrameLabel->{StyleForm["T",FontSize->18],StyleForm["<E>/n, <M>/n",FontSize->18]}]

tCV=Table[{T,CV[T]},{T,.2,5,.2}];
tchi=Table[{T,chi[T]},{T,.05,5,.2}];
plmc2=Show[ListPlot[tCV,PlotStyle->Red],ListPlot[tchi],ImageSize->600,AxesLabel->{StyleForm["T",FontSize->18],StyleForm["(Subscript[C, V])", "χ",FontSize->18]}]

Dokładne rozwiązanie Isinga

Zexa=(2Cosh[b j])^m (1+Tanh[b j]^m)

ene[T_,m_]=-D[Log[Zexa],b]/m/.j->J/.b->1/T//Simplify

CV[t_,m_]=D[ene[T,m],T]/m//FullSimplify

(* porównanie Monte Carlo i wyniku dokładnego *)
Show[plmc,Plot[ene[T,n],{T,0,5}]]
Show[plmc2,Plot[CV[T,n],{T,0,5}]]

Fluktuacje

(* ciepło właściwe porównane z teorią *)
init;
teav2=Table[{T,CV[T]},{T,.1,5,.1}];
Show[ListPlot[teav2,ImageSize->600,AxesLabel->{StyleForm["T",FontSize->18],StyleForm["(Subscript[C, V])",FontSize->18]}],Plot[CV[T,n],{T,0,5}]]

(* dokładny wynik na podatność w granicy dużych n *)
chiexa=T D[b J+ Log[Cosh[b h]+Sqrt[Sinh[b h]^2+Exp[-4b J]]],{h,2}]/.h->0 /.b->1/T //Simplify //PowerExpand

tsav2=Table[{T,chi[T]},{T,.1,5,.1}];
pl=Show[ListPlot[tsav2,PlotRange->All,ImageSize->600,AxesLabel->{StyleForm["T",FontSize->18],StyleForm["χ/n",FontSize->18]}],Plot[Exp[2J/T]/T/n,{T,.1,5},PlotRange->{0,1}]]

Korelacje

(Też można uzyskać wynik analityczny z pomocą tzw. macierzy przejścia)

gav[T_]:=Module[{},Do[sweepT[T],{k,1,ninit}];gt=Table[0,{m,0,n/3}];Do[sweepT[T];gt=gt+Table[{m,gamma[m]},{m,0,n/3}],{k,1,ncoll}];gt=gt/ncoll;ListPlot[Abs[gt],AxesLabel->{StyleForm["m",FontSize->18],StyleForm["Γ(m)",FontSize->18]},PlotRange->{-.1,1}, ImageSize->500, Joined->True]]

init;
Show[{gav[.01],gav[.8],gav[3]},PlotRange->All]

(* z konstrukcji, gamma[0]=1 *)
gamma[0]

(* dopasowanie do funkcji wykładniczej *)
gav[.8];
opt=FindFit[Abs[gt],Exp[-m g],g,m]
Show[ListPlot[Abs[gt]],Plot[Exp[-m g]/.opt,{m,0,16},AxesLabel->{StyleForm["m",FontSize->18],StyleForm["Γ(m)",FontSize->18]},PlotRange->{-.1,1}, ImageSize->500]]

ZADANIE 2

0) [onzt] Prześledź powyższy rozdział i zrób rachunek dla przypadku ferrimagnetycznego. Jakie są różnice w stosunku do ferro?
1) [onzt] Przeczytaj notatki Wei Cai (Caltech), plik Chap12_Ising_Model_v04.pdf
2) [onz2t] Rozwiąż analitycznie model Isinga (przynajmniej dla dla h=0), oddaj pisemną pracę
3) [onz3t] Do wyboru: zaprogramuj symulację Metropolisa

 - model Isinga 1D dla h≠0
 - model Isinga 1D z dwoma lub więcej rzędami
 - model Isinga 1D z oddziaływaniami dalszymi niż najbliższy sąsiad
 - model Pottsa 1D (uogólnienie Isinga dla (Subscript[s, i])=1,2,..,q), (Subscript[H, ij])=-J cos(2π((Subscript[s, i])-(Subscript[s, j]))/q)
 - model Isinga 2D
 - model Heisenberga (klasyczny) dD, (Subscript[H, ij])=-J (OverVector[(Subscript[s, i])])·(OverVector[(Subscript[s, j])]) (d wymiarów)
 
Jako wzorcem można/należy się posłużyć przykładem z tego notebooka dla modelu Isinga 1D, h=0. 
(Jeśli ktoś na prawdę woli, można uzyć inną platformę programowania, jak Matlab, Maple, itp.)

Wykład 3  
Dalsze zastosowania algorytmu Metropolisa: 
od kombinatoryki do całek po trajektoriach

Projektowanie układu scalonego

Problem jest nastepujacy: mamy n mikroelementów elektronicznych (tranzystory). Logika urządzenia określa, jaka jest sieć połączeń, tzn. które pary tranzystorów są połączone. Na jedej płytce mieści się połowa tranzystorów, więc musimy użyć dwóch płytek do realizacji układu. Połączenia tranzystrorów w obrębie płytki są tanie, natomiast połączenia między płytkami drogie (wymagają “drutów”). Podziel tranzystory między płytki w taki sposób, aby połączeń między płytkami było jak najmniej.

Rozwiązanie: jest to zakamuflowany problem spinowy. Przydziel tranzystorowi i etykietę s(i)=0 jeśli należy do pierwszej płytki, a s(i)=1 jeśli należy do drugiej. “Hamiltonian” możemy zapisać jako 

H=(∑_(i<j)) [((Subscript[s, i])+(Subscript[s, j])) mod 2].

Dalej postepujemy zgodnie z algorytmem wyżarzania, tj. po każdym przemieceniu obniżamy nieco temperaturę.

Symulacja

Quit[]

(* liczba tranzystorów *) 
n=50
(* maksymalna liczba połączeń *)
k=n(n-1)/2

ri:=Random[Integer,{1,n}]
rs:=Random[Integer,{0,1}]

(* liczba możliwości = n po n/2 *)
n!/(n/2)!/(n/2)!//N

Binomial[n,n/2] //N

(* inicjalizacja *)
init:=Module[{},t=Table[0,{n}];While[Sum[t[[i]],{i,1,n}]<n/2,t[[ri]]=1];]

init
t

(* połączenia - tutaj generujemy je losowo *)
index={};
Do[in={Sort[{ri,ri}]};index=Join[index,in],{i,1,k}]
index=DeleteDuplicates[Sort[index]];
il=Length[index]

index

energy:=Sum[Mod[t[[index[[i,1]]]]+t[[index[[i,2]]]],2],{i,1,il}]

(* operacja zmiany miejscami elementów i oraz j *)
flip[i_,j_]:=Module[{},to=t[[i]];t[[i]]=t[[j]];t[[j]]=to;]

(* przeskok Metropolisa w węźle *)
step[i_,j_,T_]:=Module[{},eold=energy;flip[i,j];enew=energy;del=enew-eold;If[del<0,nic,If[Exp[-del/T]>Random[],nic,flip[i,j]]];]

(* jedno "przemiecenie" skłąda się z n kroków *)
sweepT[T_]:=Module[{},Do[step[ri,ri,T],{i,1,n}];]

init;
t
energy
sweepT[.2]
t
energy

W naszej procedurze wyżarzania temperatura maleje z każdym przemieceniem liniowo od wartości T do wartości T/liczba przemieceń. Oczywiście, można robić granicę T ⟶ 0 w inny sposób.

etrace[T_,steps_]:=Module[{},ett={};Do[sweepT[T*(1-(k-1)/steps)];AppendTo[ett,energy],{k,1,steps}];Show[ListPlot[{ett},AxesLabel->{StyleForm["sweep",FontSize->18],StyleForm["E",FontSize->18]},ImageSize->400 ,PlotRange->All]]]

init; 
t0=t;

t=t0

t=t0
energy
etrace[.35,100]
t
energy

Wyżarzanie i szkło spinowe

(* na tablicy *)

Całki feynmanowskie po trajektoriach

Wiemy z wykładu mechaniki kwantowej, że operator ewolucji można zapisać jako całkę po trajektoriach (Feynman)

U((Subscript[t, i]),(Subscript[x, i]);(Subscript[t, f]),(Subscript[x, f])) = A ∫D[x(t)] (e^(iS((Subscript[t, i]),(Subscript[t, f]))/ℏ)),

gdzie A jest pewną stałą a S=(∫_(Subscript[t, i])^(Subscript[t, f])) dt L[x(t),(OverDot[x])(t)] jest klasycznym działaniem dla cząstki o trajektorii x(t), przy czym x((Subscript[t, i]))=(Subscript[x, i ])oraz x((Subscript[t, f]))=(Subscript[x, f ]). Przyjmujemy jednostki ℏ=1. Dyskretyzując całkę po czasie otrzymujemy, używając v=((x((Subscript[t, j+1]))-x((Subscript[t, j])))/(Δt)),

U=A ∫dx((Subscript[t, 2]))...dx((Subscript[t, n-1])) exp[iΔt(∑_(j=1)^(n-1))((((m(x((Subscript[t, j+1]))-x((Subscript[t, j]))))/(2 ((Δt)^2)))^2)-V[x((Subscript[t, j]))])],

gdzie (Subscript[t, 1])=(Subscript[t, i]), (Subscript[t, n])=(Subscript[t, f]), Δt = ((Subscript[t, f])-(Subscript[t, i]))/n. Przechodzimy teraz do urojonego czasu t=iτ:

U=A ∫(Subscript[dx, 2])...(Subscript[dx, n-1]) exp[-T(∑_(j=1)^(n-1))((((((Subscript[x, j+1])-(Subscript[x, j])))/(2m (Τ^2)/n))^2)+1/n V[(Subscript[x, j])])],

gdzie Τ=(Subscript[τ, f])-(Subscript[τ, i]), Δτ=T/n, (Subscript[x, j])=x((Subscript[t, j])). Formalnie, powyższe wyrażenie jest analogiczne do modelu spinowego w 1D, przy czym warość spinu jest zamieniona na współrzędną x i jest liczbą rzeczywistą, oddziaływanie jest między sąsiednimi węzłami (pochodzi z energii kinetycznej), a odwrotność temperatury jest proporcjonalna do euklidesowego czasu,  β∼Τ. 

Czynik normalizacji (dla przypadku swobodnego A = (((2 ((πΤ)^2)/n))^(-n/2))) nie zależy od położeń (Subscript[x, i]).

Poniżej rozważamy przypadek cząstki o masie m=1 w potencjale harmonicznym V(x).

Ogólnie, mechanika kwantowa (0 wymiarów przestrzennych,  1 wymiar czasowy, tj. zmienną jest funkcja x(t)) jest równoważna 0+1 D modelowi spinowemu, a teoria pola w d+1 wymiarach (zmienną jest pole ϕ(t,x)) jest równoważna D=d+1 wymiarowemu modelowi spinowemu (siatki QCD itd.)

PIMC - Path - Integral Monte Carlo

Quit[]

(* liczba węzłów sieci *)
n=50

(* ograniczenie na współrzędną *)
sc=3;

rn:=Random[Real,{-1,1}]sc

(* losowa inicjalizacja początkowej trajektorii *)
init:=Module[{},t=Table[rn,{i,1,n}];t[[n]]=t[[1]];]

init
t

ener[tau_]:= Sum[(t[[i]]-t[[i+1]])^2/2/(tau/n),{i,1,n-1}]/tau +Sum[t[[i]]^2/2*(tau/n),{i,1,n}]/tau

(* energia z tw. wirialnego - dla osylatora E_kin=V, E_tot=2V - taki wzor daje stabilniejsze numerycznie wyniki *)
ve[tau_]:= 2Sum[t[[i]]^2/2*(tau/n),{i,1,n}]/tau

(* operacja małej zmiany współrzędnej w węźle i *)
flip[i_]:=t[[i]]+=rn/10

(* przeskok Metropolisa w węźle i przy euklidesowym czasie tau *)
(* brutalnie zadajemy warunek t(n)=t(i) *)
step[i_,tau_]:=Module[{},tiold=t[[i]];eold=ener[tau];flip[i];enew=ener[tau];del=enew-eold;If[del<0,nic,If[Exp[-del*tau]>Random[],nic,t[[i]]=tiold]];t[[n]]=t[[1]];]

(* jedno "przemiecenie", począwszy od losowo wybranego węzła r *)
(* sweepT[tau_]:=Module[{r},r=Random[Integer,{1,n-1}];Do[step[i,tau],{i,r,n-1}];Do[step[i,tau],{i,2,r-1}];] *)

(* przemiecenie składające się z jednego kroku - większa statystyka *)
sweepT[tau_]:=Module[{r},r=Random[Integer,{1,n-1}];step[r,tau];]

Ilustracja wygładzania trajektorii

etrace[tau_,steps_]:=Module[{},ett={};Do[sweepT[tau];AppendTo[ett,ener[tau]],{k,1,steps}];Show[ListPlot[{ett},AxesLabel->{StyleForm["sweep",FontSize->18],StyleForm["E",FontSize->18]},ImageSize->400 ,PlotRange->All]]]

init;
p01=ListPlot[t,Joined->True,AxesLabel->{StyleForm["τ",FontSize->14],StyleForm["x",FontSize->14]},PlotStyle->Red,PlotRange->{-sc,sc}];
etrace[.3, 5000]
p02=ListPlot[t,Joined->True,AxesLabel->{StyleForm["τ",FontSize->14],StyleForm["x",FontSize->14]},PlotStyle->Blue,PlotRange->{-sc,sc}];

Show[p01,p02]

Symulacje

Znajdziemy dla przykładu energię i kwadrat funkcji falowej stanu  
podstawowego oscylatora harmonicznego z pomocą algorytmu Metropolisa. Oczywiście, dla oscylatora harmonicznego 
(Subscript[E, 0])=1/2 i |(Subscript[ψ, 0])((|)^2)=(e^(-(x^2)))/√π. Z drugiej strony, przechodząc do bazy energii,
 
 Z=Tr U(T)=
∫dx<x| U(Τ)|x>
=∫dx (∑_n)<x| n > < n | U(Τ) | n > < n |x>,
przy czym 
< n | U(Τ) | n > = (e^(-(Subscript[TE, n]))) 
(bo z drugiej strony w mechanice kwantowej U=(e^(-TH)), gdzie H jest niezależnym od czasu Hamiltonianem układu). 
W granicy dużych Τ mamy zatem (wyższe stany są "ubite")

∫dx< x| U |x >≃ ∫dx|(Subscript[ψ, 0])(x)((|)^2) (e^(-(Subscript[TE, 0])))=(e^(-(Subscript[TE, 0]))).

Wynika stąd, że (Subscript[E, 0])=(lim_(Τ →∞))Log[-Tr U(T)]/Τ, tj. energia swobodna dąży do energii stanu podstawowego. 

Inna metoda (zastosowana w obliczeniach poniżej):

Średnia dla niezależnego od czasu operatora O  przy ustalonym Τ wynosi 

<O(Subscript[>, Τ]) = ((Tr [O U(Τ)])/Z) = (((∑_n)<n|O|n> (e^(-Τ (Subscript[E, n]))))/((∑_n)(e^(-Τ (Subscript[E, n])))))⟶ <0|O|0> przy dużych Τ.
Ponieważ w metodzie Metropolisa stany sa próbkowane z rozkładem 
"Boltzmanna" (e^(-S)),gdzie S oznacza działanie  dla danej trajektorii,w poniżej przedstawionym programie sumy są wykonywane z wagą 1. 
A zatem w granicy dużych Τ dla O=H dostajemy energię stanu podstawowego,a dla O=|x><x| - kwadrat modułu funkcji falowej stanu podstawowego. 
Podobnie możemy liczyć inne obserwable dla stanu podstawowego.

ninit=5000; (* liczba przemieceń do osiągnięcia równowagi ... *)
ncoll=100000; (* ... i nastepująca po niej liczba przemieceń pomiarowych *)

(* średnia (po ncoll trajektoriach) energia z tw. wirialnego *)
eav2[tau_]:=Module[{},init;Do[sweepT[tau],{k,1,ninit}];scoll=0;Do[sweepT[tau];scoll+=ve[tau],{k,1,ncoll}];scoll/ncoll]

(* średnia (po ncoll trajektoriach) kwadratu modułu funkcji falowej *)
psiav2[tau_]:=Module[{},init;psi={};Do[sweepT[tau],{k,1,ninit}];Do[sweepT[tau];psi=AppendTo[psi,t],{k,1,ncoll}];psi=Flatten[psi];]

psiav2[15.]

(* symetryzacja polepsza statystykę *)
psi=Join[psi,-psi];

h=Histogram[psi,30,"PDF"]

Mean[psi]
(* wariancja powinna być równa 1/2 dla stanu podstawowego = (Subscript[E, 0 ])*)
vp=Variance[psi]

Show[h,Plot[Exp[-(x)^2]/Sqrt[Pi],{x,-4,4}]]

(* energia stanu podstawowego *)
e0=1/2

eav2[10.]

(* tablicowanie i wykres energii w funkcji τ *)
Timing[teav=Table[{tau,eav2[tau]},{tau,.5,10,.5}]]

plmc=Show[ListPlot[teav,PlotMarkers->Automatic],Plot[{e0},{tau,0,teav[[Length[teav],1]]}],ImageSize->600,Frame->True,FrameLabel->{StyleForm["τ",FontSize->20],StyleForm["(Subscript[E, 0])",FontSize->18]},PlotRange->All,AxesOrigin->{0,0}]

Mean[teav][[2]]
(StandardDeviation[teav]/Sqrt[Length[teav]])[[2]]

http://fisteo12.ific.uv.es/~santamar/qapplet/qapplet.html

Przedziały czasowe a=Τ/n muszą być mniejsze niż charakterystyczna skala, a<<1/(Subscript[E, 0]), natomiast
Τ musi być większe niż odwrotność przerwy energetycznej, Τ= n a>> 1/((Subscript[E, 1])-(Subscript[E, 0])).

Uzyskanie informacji o stanach wzbudzonych jest kosztowne.

Ulepszone działania to takie dyskretne działania, których granica ciągła a ⟶ 0 jest szybciej osiągana (tj. przy mniejszych n).

Teoria pola na siatkach: węzły - współrzędne (t, (x^1),(x^2), (x^3)), zmienne w węzłach - pola fermionowe, linki - pola cechowania.

ZADANIE 3

1) Przeczytaj skrypt “Computational Quantum Physics” de Forcranda i Troyera, rozdz. 5, 7-9.
2) Przerób i puść program z wykładu dot. całek po trajektoriach dla innego potencjału.

Wykład 4
Algorytm genetyczny i problem komiwojażera

Metoda brutalna

Quit[]

<< Combinatorica`

(* n - liczba miast do odwiedzenia, np - liczba wszystkich miast, łącznie ze startowym *)
n = 9
np = n + 1;

(* liczba wszystkich tras *)
nf = n!/2

miasto := RandomReal[{-1, 1}, 2]

(* test *)
miasto

(* tablica współrzędnych miast do odwiedzenia *)
tabmiast = 
 Table[miasto, {i, 1, n}]

(* miasto startowe *)
start = {{0, 0}}

(* tablica wszystkich miast, łącznie ze startowym *)
wszystkie = Join[start, tabmiast]

(* wykres miasta startowego *)
plstart = ListPlot[start, PlotStyle -> {PointSize[.04], RGBColor[1, 0, 1]}];

(* wykres miast do odwiedzenia *)
plotmiast = ListPlot[tabmiast, PlotStyle -> PointSize[.04]];

(* etykiety *)
lab = Table[Graphics[Text[ToString[i], tabmiast[[i]] + {.05, .05}]], {i, 1, n}];

(* nałożenie wykresów *)
pp1 = Show[plotmiast, plstart, lab, Axes -> False,PlotRange -> {{-1.1, 1.1}, {-1.1, 1.1}}]

(* odleglość euklidesowa między dwoma punktami *)
d[p1_, p2_] := Sqrt[(p1[[1]] - p2[[1]])^2 + (p1[[2]] - p2[[2]])^2]

(* tablica odległości między wszystkimi miastami *)
odleglosci =  Table[d[wszystkie[[i]], wszystkie[[j]]](* +2Random[] *), {i, 1, np}, {j, 1, np}];
odleglosci // MatrixForm

(* lista wszystkich permutacji zbioru {1,2,...,n} - sa to wszystkie możliwe \
kolejnosci odwiedzania miast *)
perm = Permutations[n];

(* długość danej drogi dr 
dr[[i]] - miasto numerowane od 0
dr[[i]]+1 - miasto numerowane od 1 *)
dlugo[dr_] := Module[{}, suma = 0; Do[ suma = suma + odleglosci[[dr[[i]] + 1]][[dr[[i + 1]] + 1]], {i, 1, np}]; suma]

Histogram[Table[dlugo[Join[{0}, perm[[Random[Integer,{1,n!}]]], {0}]],{40000}]]

(* znajdowanie najkrótszej trasy komiwojażera 
dl - długość trasy
min - aktualna najmniejsza długość trasy
i - numer trasy
imin - numer najkrótszej trasy
opti - najkrótsza trasa *)

Timing[
 Module[{}, mi = 100000; 
  Do[dl = dlugo[Join[{0}, perm[[i]], {0}]]; 
   If[dl < mi, mi = dl; imin = i], {i, 1, nf}]; 
  Print["najkrótsza droga : ", opti = perm[[imin]], ", długość : ", mi]]]

(* tablica współrzędnych optymalnej trasy, łącznie z punktem startowym i końcowym {0,0} *)

tt = Join[{{0, 0}}, Table[tabmiast[[opti[[i]]]], {i, 1, n}], {{0, 0}}];
Show[pp1, Graphics[Line[tt]], PlotRange -> All, AspectRatio -> 1]
Show[Graphics[{Purple, Opacity[.7], Polygon[tt]}], PlotRange -> All, AspectRatio -> 1]

Rozwiązanie problemu komiwojażera dla 15112 niemieckich miescowości (algorytm z powracaniem)

image = Import["/home/broniows/local/Computer_Methods/maze.gif", "GIF"]

Algorytm genetyczny

John Henry Holland

Quit[]

<< Combinatorica`

(* liczba miasteczek *)
n = 30;

(* losowe odstępstwo od euklidesowości, jesli k>0 *)
k = 0;
d[x_, y_] := Sqrt[x^2 + y^2] + k Random[Real, {0, 1}];
txy = Table[{Random[Real, {-1, 1}], Random[Real, {-1, 1}]}, {i, 1, n}];
ttu = AppendTo[txy, {0, 0}]
plot1 = ListPlot[ttu, AspectRatio -> 1, PlotStyle -> PointSize[.04], PlotRange -> {{-1.1, 1.1}, {-1.1, 1.1}}]

Dist0 =  Table[d[txy[[i]][[1]], txy[[i]][[2]]], {i, 1, n}]; (* odleglosc od (0,0) *)
Dist = Table[d[txy[[i]][[1]] - txy[[j]][[1]],txy[[i]][[2]] - txy[[j]][[2]]], {i, 1, n}, {j, 1, n}]; (* odległość dwóch punktów *)

birth[papa_, mama_] := Module[{},
  mapo = Random[Integer, {1, n - 1}]; (* miejsce podziału chromosomów *)
  papagen = Drop[papa, -(n - mapo)];
  mamagen = Drop[mama, mapo];
  If[Random[] > 0.5, mamagen = Reverse[mamagen]]; (* odwrócenie chromosomu matki *)
  mis = Complement[Union[papa], Union[papagen, mamagen]]; (* mis zawiera brakujące allele *)
  child = papagen; (* skopiuj geny ojca do dziecka *)
  Do[                                  (* korekcja chromosomu *)
   If[MemberQ[child, mamagen[[i]]], {lm = Length[mis]; (* jeśli allel się powtarza *)
     r = Random[Integer, {1, lm}]; mis[[r]]; (* to dodaj poprawiony gen mis[r] *)
     AppendTo[child, mis[[r]]]; mis = Complement[mis, {mis[[r]]}]}, 
    AppendTo[child, mamagen[[i]]]], {i, 1, Length[mamagen]}]; (* jeśli sie nie powtarza, to dodaj gen matki *)
  child]

dlugosc[path_] := 
  Dist0[[path[[1]]]] + 
   Sum[Dist[[path[[i]]]][[path[[i + 1]]]], {i, 1, n - 1}] + 
   Dist0[[path[[n]]]];

ploway[path_] := Module[{},
  tli = Table[txy[[path[[i]]]], {i, 1, n}];
  tli = AppendTo[tli, {0, 0}];
  tli = PrependTo[tli, {0, 0}];
  opti = Show[plot1, Graphics[{RGBColor[.5, .4, 0], Line[tli]}]]]

(* zapamietaj najlepszych rodziców *)
storepar[x_, wa_] := 
 Module[{}, repl = True; 
  Do[If[(x < li[[i]]) && repl, repl = False; 
    li = Drop[Insert[li, x, i], -1]; (* li - najkrótsze długości dróg rodziców *)
    bp = Drop[Insert[bp, wa, i], -1]], {i, 1, npop}]] (* bp - najlepsze drogi rodziców *)

(* zapamietaj najlepsze dzieci *)
storech[x_, wa_] := 
 Module[{}, repl = True; 
  Do[If[(x < lich[[i]]) && repl, repl = False; 
    lich = Drop[Insert[lich, x, i], -1]; (* lich - najkrótsze długości dróg dzieci *)
    bpch = Drop[Insert[bpch, wa, i], -1]], {i, 1, npop}]] (* bpch - najlepsze drogi dzieci *)

npop=20; (* liczebność populacji *)
li=Table[100000,{i,1,npop}];  (* długości dróg rodziców - inicjalizacja *)
bp=Table[{},{i,1,npop}];           (* drogi rodziców - inicjalizacja *)

ilo=7000;  (* liczebność próbki losowych dróg *)
minis=100000000;
maxis=-1;
sols={};
Timing[Do[way=RandomPermutation[n];
droga=dlugosc[way];
storepar[droga,way];
sols=Append[sols,droga];
If[droga<minis,{minis=droga;bests=way}];
If[droga>maxis,maxis=droga;worsts=way],{p,1,ilo}]]

minis (* długość najkrótszej losowej drogi *)
maxis  (* długość najdłuższej losowej drogi *)

Histogram[sols] (* rozkład długości losowych dróg *)

(*  maksymalna liczba mutacji *)
mmax = 7;
(* liczba dzieci *)
ilo = 4;
(* liczba par *)
npa = 200;
(* liczba pokoleń *)
ngen = 150;

Timing[Do[
  lich = Table[100000, {i, 1, npop}]; (* wyczyść dzieci *)
  bpch = Table[{}, {i, 1, npop}];
  Do[
   papa = bp[[Random[Integer, {1, npop}]]]; (* wybór losowy ojca *)
   mama = bp[[Random[Integer, {1, npop}]]]; (* wybór losowy matki *)
   Do[ss = Random[Integer, {-n + 1, n + 1}]; 
    child = RotateRight[ (* twórz dziecko, rodzice i dziecko są cyklicznie obracani *)
      birth[RotateLeft[papa, ss], RotateLeft[mama, ss]], ss];
    m = Random[Integer, {0, mmax}];
    Do[ (* mutacja *)
     hh = Random[Integer, {1, n}];
     child =   (* przestawienie genu z jednej losowej pozycji w inną losową pozycję *)
      Insert[Drop[child, {hh}], child[[hh]], 
       Random[Integer, {1, n}]], {u, 1, m}];
    droga = dlugosc[child]; storech[droga, child], {j, 1, ilo}], {i, 1, npa}];
nn = Floor[npop/2]; (* połowa populacji *)
  npop - nn;
  li = Join[Drop[lich, -nn], Drop[li, -(npop - nn)]];  (* uaktualnij połowę populacji rodziców dziećmi *)
  bp = Join[Drop[bpch, -nn], Drop[bp, -(npop - nn)]], {jj, 1, ngen}]]

li
ploway[bp[[1]]]  (* prawdopodobnie uwidzoczni sie degeneracja *)
dlugosc[bp[[1]]]

ZADANIE 4

1) Przeczytaj prace licencjacką W. Świderskiego (Swiderski.pdf)
2) Konkurs (nagroda dla zwycięzcy 100zł). Zaprogramuj samodzielnie w Mathematice rozwiązywanie euklidesowego problemu komiwojażera z pomocą dowolnej metody dla 200 miast o losowo wygenerowanych współrzędnych z przedziału [-1,1]. Rozstrzygnięcie konkursu nastąpi na wykładzie 9.01.2013, gdzie zostaną przez mnie podane współrzędne takie same dla wszystkich w postaci pliku w formacie {{x1,y1},{x2,y2},...,{x200,y200}}. (program musi je umieć wczytać). Każdy program będzie się jednocześnie liczył w czasie wykładu przez 60 min. i wygra ten, który znajdzie najlepsze rozwiązanie. Warunek: muszą wziąć udział co najmniej 3 programy, moga być zespołowe. Program nie może używać wbudowanych funkcji biblitecznych, np. FindShortestTour. Program może być zapożyczony, ale uczestnik konkursu musi wyjaśnić wszystkie jego kroki. Program powinien co pewien czas pokazywać aktualnie najlepsze rozwiązanie celem podgrzania emocji.

Wykład 5
Inteligencja chmary (swarm intelligence)

Algorytm kolonii mrówek (Ant Colony Optimization)

Zaproponowany przez Marco Dorigo w 1991 r. na podstawie obserwacji zachowań mrówek.

Quit[]

(* n - liczba miast do odwiedzenia *)
n = 30;

miasto := RandomReal[{-1, 1}, 2]

(* tablica współrzędnych miast do odwiedzenia *)
wszystkie = Table[miasto, {i, 1, n}];
(* wykres miasta startowego *)
plotmiast = ListPlot[wszystkie, PlotStyle -> PointSize[.04]];
(* etykiety *)
lab = Table[Graphics[Text[ToString[i], wszystkie[[i]] + {.05, .05}]], {i, 1, n}];
(* nałożenie wykresów *)
pp1 = Show[plotmiast,lab, Axes -> False,PlotRange -> {{-1.1, 1.1}, {-1.1, 1.1}},AspectRatio->1]

(* odleglość euklidesowa między dwoma punktami *)
d[p1_, p2_] := Sqrt[(p1[[1]] - p2[[1]])^2 + (p1[[2]] - p2[[2]])^2]

(* tablica odległości między wszystkimi miastami *)
odleglosci =  Table[d[wszystkie[[i]], wszystkie[[j]]], {i, 1, n}, {j, 1, n}];
odleglosci // MatrixForm

(* długość danej drogi dr *)
dlugo[dr_] := Module[{}, suma = 0; Do[suma += odleglosci[[dr[[i]]]][[dr[[i + 1]]]], {i, 1, n-1}];suma+odleglosci[[dr[[n]]]][[dr[[1]]]]]

pltr[mr_]:=Module[{},
tt=Join[Table[wszystkie[[mr[[i]]]],{i,1,n}],{wszystkie[[mr[[1]]]]}];
Show[pp1, Graphics[{Red,Line[tt]}], PlotRange -> All, AspectRatio -> 1]]

pltr0[mr_]:=Module[{},
tt=Join[Table[wszystkie[[mr[[i]]]],{i,1,n}],{wszystkie[[mr[[1]]]]}];
Show[Graphics[{Red,Line[tt]}], PlotRange -> All, AspectRatio -> 1]]

Wbudowana funkcja

Timing[mathpath=FindShortestTour[wszystkie]]

pltr0[mathpath[[2]]]

Szybciutko algorytm zachłanny

im=1;um={};
tm=odleglosci+100000IdentityMatrix[n];
Do[
km=Min[tm[[im]]];
pm=Position[tm[[im]],km][[1,1]];
um=AppendTo[um,pm];
Do[tm[[jm,im]]=100000,{jm,1,n}];im=pm,{lm,1,n}]

Print["Długość drogi zachłannej: ",dlugo[um]]
pltr0[um]

(* inicjalizacja macierzy feromonów *)
feromony=Table[1,{i,1,n},{j,1,n}]-IdentityMatrix[n];
feromony //MatrixForm

(* parametry *)
alpha=1.15;
beta=1;
h=0.03;
sc=n;
popsize=10;

nastepne[i_,doodwiedzenia_]:=Module[{p,r,s},p=Table[feromony[[i,j]]^alpha /odleglosci[[i,j]]^beta,{j,doodwiedzenia}];p=p/Apply[Plus,p];s=Accumulate[p];r=Random[];jg=1;While[r>s[[jg]],jg++];doodwiedzenia[[jg]]]

one[pop_]:=Module[{},fer2=Table[0,{i,1,n},{j,1,n}];dl0=Infinity;
Do[
is=Random[Integer,{1,n}]; (* miasto startowe dla danej mrowki *)
mr={is};
dood=Join[Table[i,{i,1,is-1}],Table[i,{i,is+1,n}]];

Do[jn=nastepne[is,dood];is=jn;
AppendTo[mr,jn];
pd=Position[dood,jn][[1,1]];
dood=Drop[dood,{pd,pd}],{k,2,n}];

dl=dlugo[mr]; If[dl<dl0,dl0=dl;best0=mr];

Do[fer2[[mr[[k]],mr[[k+1]]]]+= sc/dl,{k,1,n-1}];fer2[[mr[[n]],mr[[1]]]]+= sc/dl,{st,1,pop}]]

dl00=Infinity;

Timing[Do[one[popsize];If[dl0<dl00,dl00=dl0;best00=best0];
fer2=fer2/popsize;
feromony=(1-h)feromony+h*fer2;
feromony=(feromony+Transpose[feromony])/2,{st,1,500}]]

(* Chop[feromony]//MatrixForm *)
Print["najkrótsza znaleziona do tej pory: ",dl00]
pltr0[best00]

Show[plotmiast,Table[Graphics[{AbsoluteThickness[7feromony[[i,j]]],Line[{wszystkie[[i]],wszystkie[[j]]}]}],{i,1,n},{j,1,n}],pltr[best00],Axes->False,AspectRatio->1,ImageSize->500]

ZADANIE 5

1) Pobaw się parametrami algorytmu mrówkowego i opisz jakościowo zachowanie w zależności od parametrów α i β

Wykład 6
Sieci neuronowe I

Zamiast wstępu

Złudzenie Macha

Quit[]

(* funkcja intensywności *)
low=0.2; up=0.8;
xl=300; xu=600;

p[x_]:=low /; x< xl
p[x_]:=up/; x>xu
p[x_]:= low+(x-xl)/(xu-xl)*(up-low) /;x ≥xl && x≤ xu

Plot[p[x],{x,0,900},PlotRange->{0,1},Frame->True,FrameLabel->{StyleForm["x",FontSize->18],StyleForm["intenstywność",FontSize->18]}]

pic=Transpose[Table[{p[i]},{i,0,900}]]

ArrayPlot[pic,PlotRange->{0,1},ImageSize->900,AspectRatio->1,Frame->False]

Pojawienie się jaśniejszego i ciemniejszego paska jest złudzeniem, bo funkcja intensywności jest monotoniczna. Dowodzi to “korekcji” obrazu, zachodzącej w mózgu.

Złudzenie Hermanna

width=5;gap=1;nsquares=6;
Show[Graphics[hermann=Flatten[Table[{Rectangle[{x,y},{x+width,y+width}]},{x,0,(width+gap)*(nsquares-1),width+gap},{y,0,(width+gap)*(nsquares-1),(width+gap)}],1]],AspectRatio->1,ImageSize->900]

Fala liści

Show[Import["/home/broniows/kielce/Computer_Methods_2015/iluzje/leafwave.jpg","JPG"],ImageSize->1400]

Rotujące węże

Show[Import["/home/broniows/kielce/Computer_Methods_2015/iluzje/rotsnake.gif","GIF"],ImageSize->1300]

Spostrzeżenia

System wzrokowy (siatkówka, mózg) dokonuje korekcji obrazu
Efekt jest dynamiczny (falowanie), tj. korygowany na bieżąco (NB skala czasowa jest niespodziewanie długa, orzędu 1s)
Ponieważ zachodzi dla każdego, jest na poziomie głebokim (wetware)

Sieci neuronowe - podstawowe informacje

Show[GraphicsGrid[{{Show[Import["/home/broniows/local/Computer_Methods/neuron/neuron.png","PNG"],ImageSize->500],
Show[Import["/home/broniows/local/Computer_Methods/neuron/neural_network.png","PNG"],ImageSize->300]}}]]

Najprostsza sieć “w przód”: warstwa inputu (dendryty), ukryte warstwy pośrednie (może być ich wiele), warstwa outputu (aksony). Połączenia między sąsiednimi warstwami (akson-dendryt) - jest ich bardzo wiele. 

Mózg człowieka zawiera ok. (10^11) neuronów, kazdy ma średnio 7000 połączeń z innymi neuronami. Pojedynczy neuron działa niezwykle prosto (przekaźnik progowy), jednak zachowanie całego mózgu jest bardzo złożone. 

“Złożoność z prostoty”.

Prosty model neuronu (Warren McCulloch, Walter Pitts, 1942)

Show[Import["/home/broniows/kielce/Computer_Methods_2015/neuron/McCulloch-Pitts.png","PNG"],ImageSize->600]

(Subscript[x, i]) - sygnały wejścia, (Subscript[w, i]) ∈ℝ - wagi, (Subscript[w, 0 ])- uprzedzenie (ang. bias), f - funkcja progowa
(Subscript[w, i])>0 - wzbudzenie, (Subscript[w, i])<0 - tłumienie

f[S_]=UnitStep[S];
Plot[f[S],{S,-2,2},PlotStyle->Thickness[.01],AxesLabel->{StyleForm["sygnał",FontSize->18],StyleForm["nie odpali / odpali",FontSize->18]},ImageSize->600]

Zauważmy, że funkcja f jest nieliniowa, tzn. f(x+y) ≠ f(x)+f(y). Nieliniowość jest kluczową własnością układów neuronowych. Ponieważ nieciągłość nie jest fizyczna, często dokonuje się wygładzenia funkcji wzgidzenia, stosując tzw. “sigmoid”:

fs[S_]=1/(1+Exp[-S])

Plot[{fs[S],Tanh[S/2]/2+1/2},{S,-8,8},PlotStyle->Thickness[.01],AxesLabel->{StyleForm["S",FontSize->18],StyleForm["sigmoid=tanh(S/2)+1/2",FontSize->18]},ImageSize->600]

Model neuronu MCP (z ostrym progiem):

neuron[{x1_,x2_},{w0_,w1_,w2_}]=f[w0+w1 x1+w2 x2]

neuron[{.2,-.3},{-1,2,1}]
neuron[{3,.5},{-1,2,1}]

Bramki logiczne i obwody logiczne

Z pomocą neuronu MCP możemy tworzyć sieci logiczne.

(* bramka AND *)
neuron[{0,0},{-1,.6,.6}]
neuron[{1,0},{-1,.6,.6}]
neuron[{0,1},{-1,.6,.6}]
neuron[{1,1},{-1,.6,.6}]

neurAND[x1_,x2_]=neuron[{x1,x2},{-1,.6,.6}];

(* bramka NAND *)
neuron[{0,0},{1,-0.6,-0.6}]
neuron[{1,0},{1,-0.6,-0.6}]
neuron[{0,1},{1,-0.6,-0.6}]
neuron[{1,1},{1,-0.6,-0.6}]

neurNAND[x1_,x2_]=neuron[{x1,x2},{1,-0.6,-0.6}];

(* bramka OR *)
neuron[{0,0},{-1,1.2,1.2}]
neuron[{1,0},{-1,1.2,1.2}]
neuron[{0,1},{-1,1.2,1.2}]
neuron[{1,1},{-1,1.2,1.2}]

neurOR[x1_,x2_]=neuron[{x1,x2},{-1,1.2,1.2}];

Rotate[Show[Import["/home/broniows/kielce/Computer_Methods_2015/neuron/xor.jpg","JPG"],ImageSize->400],Pi/2]

(* idziemy w tył do kolejnych warstw sieci *)
neurXOR[x1_,x2_]=neurAND[neurNAND[x1,x2],neurOR[x1,x2]];

(* bramka XOR *)
neurXOR[0,0]
neurXOR[1,0]
neurXOR[0,1]
neurXOR[1,1]

neurNOT[x_]=neuron[{x,0},{-1.2x +0.2,0.2,0}]

neurNOT[0]
neurNOT[1]

Widzimy więc, że jesteśmy w stanie utworzyć podstawowe bramki logiczne. Z bramek tych możemy z kolei utworzyć dowolną sieć logiczną (bramki NAND oraz bramki NOR są zupełne, tj. każdą sieć można zrealizować w oparciu jedynie o bramki NAND lub bramki NOR). Tak więc wielowarstwowa sieć z bramek MCP jest równoważna uniwersalnej maszynie liczącej Turinga.

Rozpoznawanie wzorców

Show[Import["/home/broniows/kielce/Computer_Methods_2015/neuron/TH.png","PNG"],ImageSize->1200]

Reguła Hamminga - najkrótsza "odległość" od wzorca

górny wiersz                środkowy wiersz               dolny wiersz

X11	1  1  1  0  1  0  0  0      X21    0  1    1   1   0    1    0    0      X31   0  1   1  1  0  1  0  0
X12	1  0  1  1  0  1  0  0      X22    1  1    1   0   1    0    0    0      X32   1  0   1  1  1  0  0  0
X13	1  1  0  1  0  0  1  0      X23    0  1    0   1   1    0    1    0      X33   0  1   1  0  1  0  1  0

F1	1  0  1  1  0  1  0  0      F2      1  0   0/1 0   0/1 0/1 0/1 1      F3     1  0   0  1  1  0  0  1

Kolumny przy Xij odpowiadają wierszom na rysunku, pierwsze kolumy są dla T, drugie dla H. Wartości przy Fi są odpowiednimi wartościami funkcji Hamminga.

aa={1,2,3}

aa//FullForm

aa[[0]]

ff=Apply[g,aa]//FullForm

d[s1_,s2_]:=Apply[Plus,(s1-s2)^2]  (* odległość łańcuchów *)

f1[s_]:=Which[d[{1,1,1},s]<d[{1,0,1},s],1,d[{1,1,1},s]>d[{1,0,1},s],0,d[{1,1,1},s]==d[{1,0,1},s],1/2]
f2[s_]:=Which[d[{0,1,0},s]<d[{1,1,1},s],1,d[{0,1,0},s]>d[{1,1,1},s],0,d[{0,1,0},s]==d[{1,1,1},s],1/2]
f3[s_]:=Which[d[{0,1,0},s]<d[{1,0,1},s],1,d[{0,1,0},s]>d[{1,0,1},s],0,d[{0,1,0},s]==d[{1,0,1},s],1/2]

{Transpose[Flatten[Table[{i,j,k,f1[{i,j,k}]},{i,0,1},{j,0,1},{k,0,1}] ,2]] //MatrixForm,
Transpose[Flatten[Table[{i,j,k,f2[{i,j,k}]},{i,0,1},{j,0,1},{k,0,1}] ,2]] //MatrixForm,
Transpose[Flatten[Table[{i,j,k,f3[{i,j,k}]},{i,0,1},{j,0,1},{k,0,1}] ,2]] //MatrixForm}

g[s1_,s2_,s3_]:=Graphics[Raster[{0.9-0.4s3,0.9-0.4s2,0.9-0.4s1}]]

Dzielimy odpowiedzi na klasy w zależności od odległości od wzorca:

g[{1,1,1},{0,1,0},{0,1,0}]

rozpoznaj[s1_,s2_,s3_]:=Module[{},test=f1[s1]+f2[s2]+f3[s3];Print[GraphicsGrid[{{g[s1,s2,s3],g[{f1[s1]},{f2[s2]},{f3[s3]}]}}],Which[test==3,"to jest T",test<3 && test >=2,"raczej T", test>1 && test<2,"nie wiem",test≤1 && test >0,"raczej H",test==0,"to jest H"]]]

rozpoznaj[{1,1,1},{0,1,0},{0,1,0}]
rozpoznaj[{1,0,1},{1,1,1},{1,0,1}]
rozpoznaj[{1,0,1},{1,0,0},{1,0,1}]

r:=Random[Integer,{0,1}]

xr1={r,r,r}
xr2={r,r,r}
xr3={r,r,r}
g[xr1,xr2,xr3]

rozpoznaj[xr1,xr2,xr3]

Do[rozpoznaj[{r,r,r},{r,r,r},{r,r,r}],{20}]

Neuron liniowy jako model pamięci

Quit[]

A=
       {{0,0,0,0,0,0,0,0,0,0,0,0},
     {0,0,0,0,0,1,0,0,0,0,0,0},
     {0,0,0,0,0,1,0,0,0,0,0,0},
     {0,0,0,0,1,1,1,0,0,0,0,0},
     {0,0,0,0,1,0,1,0,0,0,0,0},
     {0,0,0,1,0,0,0,1,0,0,0,0},
     {0,0,0,1,0,0,0,1,0,0,0,0},
     {0,0,1,1,1,1,1,1,1,0,0,0},     
     {0,0,1,0,0,0,0,0,1,0,0,0},
     {0,1,0,0,0,0,0,0,0,1,0,0},  
     {0,1,0,0,0,0,0,0,0,1,0,0},
     {0,0,0,0,0,0,0,0,0,0,0,0}};

a=
       {{0,0,0,0,0,0,0,0,0,0,0,0},
     {0,0,0,0,0,0,0,0,0,0,0,0},
     {0,0,0,0,0,0,0,0,0,0,0,0},
     {0,0,0,0,1,1,1,1,0,0,0,0},
     {0,0,0,1,0,0,0,1,0,0,0,0},
     {0,0,0,1,0,0,0,1,0,0,0,0},
     {0,0,0,1,0,0,0,1,0,0,0,0},
     {0,0,0,0,1,1,1,1,1,0,0,0},     
     {0,0,0,0,0,0,0,0,0,0,0,0},
     {0,0,0,0,0,0,0,0,0,0,0,0},  
     {0,0,0,0,0,0,0,0,0,0,0,0},
     {0,0,0,0,0,0,0,0,0,0,0,0}};

ii=
       {{0,0,0,0,0,0,0,0,0,0,0,0},
     {0,0,0,0,0,0,0,0,0,0,0,0},
     {0,0,0,0,0,0,0,0,0,0,0,0},
     {0,0,0,0,0,0,1,0,0,0,0,0},
     {0,0,0,0,0,0,0,0,0,0,0,0},
     {0,0,0,0,0,1,1,0,0,0,0,0},
     {0,0,0,0,0,0,1,0,0,0,0,0},
     {0,0,0,0,0,0,1,0,0,0,0,0},     
     {0,0,0,0,0,0,1,0,0,0,0,0},
     {0,0,0,0,0,0,0,1,1,0,0,0},  
     {0,0,0,0,0,0,0,0,0,0,0,0},
     {0,0,0,0,0,0,0,0,0,0,0,0}};

II=
      {{0,0,0,0,0,0,0,0,0,0,0,0},
     {0,0,0,0,1,1,1,0,0,0,0,0},
     {0,0,0,0,0,1,0,0,0,0,0,0},
     {0,0,0,0,0,1,0,0,0,0,0,0},
     {0,0,0,0,0,1,0,0,0,0,0,0},
     {0,0,0,0,0,1,0,0,0,0,0,0},
     {0,0,0,0,0,1,0,0,0,0,0,0},
     {0,0,0,0,0,1,0,0,0,0,0,0},     
     {0,0,0,0,0,1,0,0,0,0,0,0},
     {0,0,0,0,0,1,0,0,0,0,0,0},  
     {0,0,0,0,1,1,1,0,0,0,0,0},
     {0,0,0,0,0,0,0,0,0,0,0,0}};

cosik=
      {{0,1,0,0,0,0,0,0,0,0,0,0},
     {0,0,1,0,0,0,0,0,0,0,0,0},
     {0,0,0,1,0,0,0,0,0,0,1,0},
     {0,0,0,0,0,0,0,0,0,1,0,0},
     {0,0,0,0,0,0,0,0,1,0,0,0},
     {0,0,0,0,0,0,0,1,0,0,0,0},
     {0,0,0,0,0,0,1,0,0,0,0,0},
     {0,0,0,0,0,1,0,0,0,0,0,0},     
     {0,0,0,0,0,1,0,0,0,0,0,0},
     {0,0,0,0,0,1,0,0,0,0,0,0},  
     {0,0,0,0,1,1,1,1,1,1,0,0},
     {0,0,0,0,0,0,0,0,0,0,0,0}};

orig=GraphicsGrid[{{ArrayPlot[A],ArrayPlot[a],ArrayPlot[II],ArrayPlot[ii],ArrayPlot[cosik]}},ImageSize->800]

Rozwijamy powyższe tablice w wektory celem łatwiejszej obróbki algebraicznej:

fA=Flatten[A];
fa=Flatten[a];
fI=Flatten[II];
fi=Flatten[ii];
fcosik=Flatten[cosik];

Pamięć heteroasocjatywna (kojarzeniowa)

Dygresja z algebry liniowej:
Iloczyn zewnetrzny (lub diadyczny) wektorów tworzy macierz (diadę), np.

v={v1,v2,v3};
w={w1,w2};
m=Outer[Times,v,w];
m // MatrixForm

m.w /w.w //Simplify //MatrixForm

vp={vp1,vp2,vp3};
wp={wp1,wp2};
m=Outer[Times,v,w]+Outer[Times,vp,wp];
m // MatrixForm

m.w/w.w-v //Simplify //MatrixForm

m.w / w.w   //FullSimplify //MatrixForm
m.w / w.w  -v //FullSimplify //MatrixForm
m.w / w.w  -(v +vp w.wp/w.w) //Simplify

m.wp / wp.wp  //FullSimplify //MatrixForm
m.wp / wp.wp  -vp //FullSimplify //MatrixForm
m.wp / wp.wp  -(vp +v w.wp/wp.wp) //Simplify

Uzyskujemy porządane własności 

m w / w.w = v, 
m w’ / w’.w’ = v’, 

jeśli w.w’=0, tzn. jeśli wektory w i w’ są ortogonalne, tj.

cos α = w.w’ / (Sqrt[w.w w'.w']) = 0.

Dla dużej liczby wymiarów n przypadkowo wybrane k wektory, k<<n, bedą w przybliżeniu ortogonalne. Ta własność jest użyta w konstrukcji rozpoznawania wzorców ukazanej poniżej.

fa
fI
fa.fI /Sqrt[fa.fa fI.fI]//N

Tworzymy teraz “liniowy model pamięci skojarzeniowej” (heteroasocjacja.
Partition działa odwrotnie do Flatten i z wektora robi z powrotem tablice.

W=Outer[Times,fa,fA]+Outer[Times,fA,fa]+Outer[Times,fi,fI]+Outer[Times,fI,fi];
pl[x_]:=ArrayPlot[Partition[W.x,12]/x.x,PlotRange->{0,2.5}]

Uzyskujemy skojarzenia między parami a-A, i-I. “Przekłamania” wynikają z nieortogonalności

orig
GraphicsGrid[{{pl[fA],pl[fa],pl[fI],pl[fi],pl[fcosik]}},ImageSize->800]

Spróbujmy pozbyć się przekłamań poprzez założenie odpowiedniego filtra, uzuwającego “najjaśniejsze odcienie”. Floor[x] daje największą liczbe całkowitą mniejszą lub równą x. Dodatkowo, wykonujemy Unitize, co wszystkie liczby całkowite >0 zastepuje liczbą 1. Np. dla “i” mamy

Partition[W.fi,12]/fi.fi //N //MatrixForm
Floor[%] //MatrixForm
Unitize[%] //MatrixForm

Założenie naszego (nieliniowego!) filtra w ewidently sposób pomaga! Stworzylismy zatem model pamieci heteroasocjatywnej.

orig
plf[x_]:=ArrayPlot[Partition[W.x,12]/x.x //Floor//Unitize]
GraphicsGrid[{{plf[fA],plf[fa],plf[fI],plf[fi],plf[fcosik]}},ImageSize->800]

Pamięć autoasocjatywna

Teraz kojarzymy obrazy same ze sobą (autoasocjacja)

W=Outer[Times,fa,fa]+Outer[Times,fA,fA]+Outer[Times,fi,fi]+Outer[Times,fI,fI];

pl[x_]:=ArrayPlot[Partition[W.x,12]/x.x]
plf[x_]:=ArrayPlot[Partition[W.x,12]/x.x //Floor //Unitize]

Niezaszumiony obraz

Show[GraphicsGrid[{{pl[fi],plf[fi]}}],ImageSize->400]
Show[GraphicsGrid[{{pl[fA],plf[fA]}}],ImageSize->400]
Show[GraphicsGrid[{{pl[fI],plf[fI]}}],ImageSize->400]
Show[GraphicsGrid[{{pl[fa],plf[fa]}}],ImageSize->400]

Zaszumiony obraz, tj. liczba “zle” bitów jest przypadkowo dolosowana w losowych miejscach. Rozpoznawanie działa dobrze, jeśli przekłamań nie jest zbyt wiele. Istotny jest próg filtra.

zle=10;
zn=A;
fAcor=Flatten[zn];
Do[poz=Max[fAcor];i=Random[Integer,{1,Length[fAcor]}];fAcor[[i]]=Random[Integer,{0,1}],{zle}]
GraphicsGrid[{{ArrayPlot[Partition[fAcor,12]],pl[fAcor],plf[.88 poz fAcor],plf[.5 poz fAcor]}},ImageSize->600]

cor[s_,zl_,po_]:=Module[{fAcor},fAcor=Flatten[s];
Do[i=Random[Integer,{1,Length[fAcor]}];fAcor[[i]]=Random[Integer,{0,1}],{zl}]
GraphicsGrid[{{ArrayPlot[Partition[fAcor,12]],pl[fAcor],plf[po fAcor]}},ImageSize->600]]

cor[A,5,.7]
cor[a,5,.7]
cor[II,5,.7]
cor[ii,5,.7]

Uzyskujemy korekcję błędu z zaszumionego/niekompletnego inputu

zgadnij[zl_]:=Module[{},si=Random[Integer,{1,4}];Which[si==1,s=A,si==2,s=a,si==3,s=II,si==4,s=ii];fAcor=Flatten[s];
Do[i=Random[Integer,{1,Length[fAcor]}];fAcor[[i]]=Random[Integer,{0,1}],{zl}];
GraphicsGrid[{{ArrayPlot[Partition[fAcor,12]] (*,pl[fAcor],plf[po fAcor]*)}},ImageSize->200]]

zgadnij[70]

GraphicsGrid[{{ArrayPlot[Partition[fAcor,12]] ,pl[fAcor] (*,plf[po fAcor]*)}},ImageSize->400]

p=.66;
GraphicsGrid[{{ArrayPlot[Partition[fAcor,12]] ,pl[fAcor],plf[p fAcor]}},ImageSize->600]

Which[si==1,"A",si==2,"a",si==3,"I",si==4,"i"]

Dygresja o sieciach liniowych

Zobaczyliśmy właśnie, że sieci liniowe (w szczególności po założeniu na końcu obliczeń nieliniowego filtra) potrafią wykonywać użyteczne obliczenia. Struktura tych sieci w ujeciu neuronowym wygląda następująco: mają dwie warstwy, input i output, których wierzchołki połączone są krawędziami z wagami (Subscript[w, ij]). Działanie jest tożsame z mnożeniem macierzowym

(Subscript[y, i]) = (∑_j)(Subscript[w, ij])(Subscript[x, j]).

Ze względu na liniowość nie ma sensu rozważać sieci liniowych o większej liczbie warstw, bowiem wówczas mielibyśmy

(Subscript[y, i]) = (∑_k)(Subscript[(w^((1))), ik])(∑_j)(Subscript[(w^((2))), kj])(Subscript[x, j]) = (∑_j)(Subscript[w, ij])(Subscript[x, j]),

gdzie (Subscript[w, ij]) = (∑_k)(Subscript[(w^((1))), ik])(Subscript[(w^((2))), kj]) . 

Sieci liniowe wielowarstwowe (z ew. założonym na ostatniej warswie filtrem, np. “sigmoid”) są zatem równoważne sieci dwuwarstwowej input-output. Sieci te nie mogą wykonywać wszystkich obliczeń, w szczególności nie można z ich pomocą zaimplementować bramki XOR.

Część kory mózgowej odpowiedzialna za postrzeganie wzrokowe

Show[Import["/home/broniows/kielce/Computer_Methods_2015/neuron/Brodmann.png","PNG"],ImageSize->500]

HCP - Human Connectome Project http://www.humanconnectomeproject.org/

ZADANIE 6

1) Przeczytaj artykuł  Martijn P. van den Heuvel, Olaf Sporns, “Rich-Club Organization of the Human Connectome”, The Journal of Neuroscience, 31(44):15775 (2011) - plik 15775.full.pdf

Wykład 7
Sieci neuronowe II

Uczenie nadzorowane sieci neuronowej

Perceptron (Frank Rosenblatt, 1957)

warstwa inputowa                          warstwy ukryte                           warstwa outputowa

Najprostsza sieć neuronowa z warstwą inputu i pojedynczą (lub wielokrotną) warstwą neuronów.
NB, opublikowanie książki “Perceptrons: an introduction to computational geometry” (Marvin Minsky, Seymour Papert, 1969) doprowadziło do “zimy sztucznej inteligencji”. Sieci wielowarstwowe nie mają omawianych tam ograniczeń.

Dwa etapy: 1) uczenie 2) używanie

Niech sieć o wagach {w} przeprowadza wektor inputowy x w wektor outputowy y, tj. y=S({w},x)
Niech pary (x,t) oznaczają input i oczekiwany (właściwy) output (musimy znać odpowiedzi wcześniej). Par tych jest N i stanowią one naszą próbkę szkoleniową.
  
1) Dla danej próbki szkoleniowej przeprowadzamy uczenie, polegające na takim doborze wag, aby błąd był jak najmniejszy. Robimy to np. tak: w pętli a=1..N bierzemy (Subscript[x, a]) i propagujemy wprzód, otrzymując (Subscript[y, a])=S({w},(Subscript[x, a])). Obliczamy błąd, zazwyczaj postaci 

E({w})=(∑_a)((((Subscript[y, a]) - (Subscript[t, a])))^2).

Błąd ten zależy oczywiście od doboru wag. Traktujemy teraz błąd jako funkcję {w} (dla ustalonej próbki szkoleniowej) i po prostu minimalizujemy E({w}). Nie jest to jednak procedura łatwa z powodu wystepowania minimów lokalnych. 

Popularnym algorytmem jest przedstawiona poniżej, oparta o metodę najszybszego spadku, propagacja wsteczna, polegająca na uaktualnianiu wag idąc od tyłu, skąd nazwa.

2) Ustaliwszy wagi, mamy sieć dedykowaną danemu zadaniu. Dla danego x uzyskujemy “odpowiedź” y=S({(Subscript[w, min])},x). 

Uwagi: 
Perceptron jest pewną nieliniową interpolacją. 
“Programowanie” polega na znalezieniu odpowiednich wag.
Lokalne minima są istotnym problemem - naiwny algorytm uczenia nie zawsze zadziała.
Analogia do innych sieci, skąd można brać bardziej wyszukane metody minimalizacji (wyżarzanie, algorytm genetyczny, mrówkowy,...)
Sieci o wiekszej liczbie warstw mają wieksze możliwości “rozpoznawania”, ale wymagaja dłuższego uczenia i są bardziej podatne na problem lokalnych minimów.

Sigmoid

Quit[]

(* sigmoid *)
f[x_]=1/(1+Exp[-x])  (* rozkład Fermiego-Diraca z odwrotnym znakiem x :) *)

Plot[f[x],{x,-6,6}]

a1=D[f[x],x] 
a2=f[x](1-f[x])//Simplify

a1/a2 //Simplify

Dla sigmoidu f’(x)=f(1-f)

a1==a2 //Simplify

Niekiedy wprowadza się sigmoid z “temperaturą”, f(x;T)=1/(1+Exp[-x/T]). Parametr T rozmywa “obszar decyzyjny”. W naszych zastosowaniach T=1.

Propagacja wsteczna

Wykład 7
Sieci neuronowe II

Uczenie nadzorowane sieci neuronowej

Perceptron (Frank Rosenblatt, 1957)

warstwa inputowa                          warstwy ukryte                           warstwa outputowa

Najprostsza sieć neuronowa z warstwą inputu i pojedynczą (lub wielokrotną) warstwą neuronów.
NB, opublikowanie książki “Perceptrons: an introduction to computational geometry” (Marvin Minsky, Seymour Papert, 1969) doprowadziło do “zimy sztucznej inteligencji”. Sieci wielowarstwowe nie mają omawianych tam ograniczeń.

Dwa etapy: 1) uczenie 2) używanie

Niech sieć o wagach {w} przeprowadza wektor inputowy x w wektor outputowy y, tj. y=S({w},x)
Niech pary (x,t) oznaczają input i oczekiwany (właściwy) output (musimy znać odpowiedzi wcześniej). Par tych jest N i stanowią one naszą próbkę szkoleniową.
  
1) Dla danej próbki szkoleniowej przeprowadzamy uczenie, polegające na takim doborze wag, aby błąd był jak najmniejszy. Robimy to np. tak: w pętli a=1..N bierzemy (Subscript[x, a]) i propagujemy wprzód, otrzymując (Subscript[y, a])=S({w},(Subscript[x, a])). Obliczamy błąd, zazwyczaj postaci 

E({w})=(∑_a)((((Subscript[y, a]) - (Subscript[t, a])))^2).

Błąd ten zależy oczywiście od doboru wag. Traktujemy teraz błąd jako funkcję {w} (dla ustalonej próbki szkoleniowej) i po prostu minimalizujemy E({w}). Nie jest to jednak procedura łatwa z powodu wystepowania minimów lokalnych. 

Popularnym algorytmem jest przedstawiona poniżej, oparta o metodę najszybszego spadku, propagacja wsteczna, polegająca na uaktualnianiu wag idąc od tyłu, skąd nazwa.

2) Ustaliwszy wagi, mamy sieć dedykowaną danemu zadaniu. Dla danego x uzyskujemy “odpowiedź” y=S({(Subscript[w, min])},x). 

Uwagi: 
Perceptron jest pewną nieliniową interpolacją. 
“Programowanie” polega na znalezieniu odpowiednich wag.
Lokalne minima są istotnym problemem - naiwny algorytm uczenia nie zawsze zadziała.
Analogia do innych sieci, skąd można brać bardziej wyszukane metody minimalizacji (wyżarzanie, algorytm genetyczny, mrówkowy,...)
Sieci o wiekszej liczbie warstw mają wieksze możliwości “rozpoznawania”, ale wymagaja dłuższego uczenia i są bardziej podatne na problem lokalnych minimów.

Sigmoid

Quit[]

(* sigmoid *)
f[x_]=1/(1+Exp[-x])  (* rozkład Fermiego-Diraca z odwrotnym znakiem x :) *)

Plot[f[x],{x,-6,6}]

a1=D[f[x],x] 
a2=f[x](1-f[x])//Simplify

a1/a2 //Simplify

Dla sigmoidu f’(x)=f(1-f)

a1==a2 //Simplify

Niekiedy wprowadza się sigmoid z “temperaturą”, f(x;T)=1/(1+Exp[-x/T]). Parametr T rozmywa “obszar decyzyjny”. W naszych zastosowaniach T=1.

Propagacja wsteczna

http://www.speech.sri.com/people/anand/771/html/node37.html

Skonstruujemy teraz perceptron o 2 warstwach ukrytych

ni=3; (* liczba neuronów wejściowych *)
n1=7;  (* liczba neuronów pierwszej warstwy ukrytej *)
n2=7;  (* liczba neuronów drugiej warstwy ukrytej *)
no=1;   (* jeden neuron wyjściowy *)

Warunki dla trzech rozważanych przykładów - trzeba wczytać tylko jeden (!)

Cond[x_]:=x[[1]]^2+x[[2]]^2<0.75^2 
cas=1;

Cond[x_]:=x[[1]]>x[[2]]
cas=2;

Cond[x_]:=(x[[1]]-.5)^2+(x[[2]]-.5)^2<0.45^2
cas=3;

(* ten warunek "wpada" w lokalne minimum i mój układ nie działa! *)
Cond[x_]:=(x[[1]]-.5)^2+(x[[2]]-.5)^2<0.45^2 && (x[[1]]-.5)^2+(x[[2]]-.5)^2>0.3^2
cas=4;

cas

Początkowe wagi wybrane sa losowo z przedziału (-et1, et1)

et1=2.25;

W1=Table[et1(1-2Random[]),{ni},{n1}]
W2=Table[et1(1-2Random[]),{n1},{n2}]
W3=Table[et1(1-2Random[]),{n2},{no}]

Algorytm wstecznej propagacji błędu (Arthur E. Bryson, Yu-Chi Ho, 1969. Od połowy lat 70. został wprowadzony z powodzeniem do sieci neuronowych i spowodował renesans badań nad sztuczną inteligencją.

E=1/2(((t-y))^2)
y = f(∑(Subscript[x, i])(Subscript[w, i]))

((∂f(s))/(∂s))=f(s)(1-f(s))                   - sigmoid

((∂E)/(∂(Subscript[w, i])))=((∂E)/(∂y))((∂y)/(∂s))((∂s)/(∂(Subscript[w, i])))=(y-t)y(1-y)(Subscript[x, i])

nstep=50000; (* liczba kroków uczenia musi być dostatecznie duża *)

Do[x={Random[],Random[],1};
If[Cond[x],t=1,t=0]; (* output ma dać 1 lub 0 *)

eta=0.02; (* prędkość uczenia *)

y1=f[x.W1]; (* wektory przepływu w kolejnych warstwach *)
y2=f[y1.W2];
y3=f[y2.W3];

(* propagacja wsteczna wg. algorytmu Bryson-Ho *)

D3=Table[y3[[i]](1-y3[[i]])(t-y3[[i]]),{i,1,no}];
W3+=Table[eta D3[[i]] y2[[j]],{j,1,n2},{i,1,no}];
D2=Table[y2[[i]](1-y2[[i]])Sum[W3[[i,j]]D3[[j]],{j,1,no}],{i,1,n2}];
  W2+=Table[eta D2[[i]] y1[[j]],{j,1,n1},{i,1,n2}];
D1=Table[y1[[i]](1-y1[[i]])Sum[W2[[i,j]]D2[[j]],{j,1,n2}],{i,1,n1}];
  W1+=Table[eta D1[[i]]x[[j]],{j,1,ni},{i,1,n1}],{nstep}]

Po nauce czas na użycie ...

(* pojedyncze uzycie *)
x={Random[],Random[],1}
Cond[x]
y1=f[x.W1];
y2=f[y1.W2];
y3=f[y2.W3][[1]]

(Tu powrócić później)

(* przecięcie pewnych neuronów *)
W1[[1,2]]=0;
W1[[1,3]]=0;
W1[[1,4]]=0;

sampP:=Module[{},x={Random[],Random[],1};
y1=f[x.W1];
y2=f[y1.W2];
y3=f[y2.W3][[1]];{If[y3>0.5,Red,Blue],Point[{x[[1]],x[[2]]}]}]

Switch[cas,
1,gr=Plot[Sqrt[.75^2-x^2],{x,0,1},PlotStyle->Thickness[.004]],
2,gr=Plot[x,{x,0,1},PlotStyle->Thickness[.004]],
3,gr=ParametricPlot[{.45 Cos[th]+.5,.45Sin[th]+.5},{th,0,2Pi},PlotStyle->Thickness[.004]]];

Show[Table[Graphics[sampP],{2000}],gr,ImageSize->700]

Quiz

A jak my sobie radzimy “na oko”?

x0={Random[],Random[],1}

x0={0.367126,0.935533,1}

Show[Plot[-5,{x,0,1},PlotRange->{{0,1},{0,1}}],Graphics[{PointSize->.01,Point[{x0[[1]],x0[[2]]}]}],ImageSize->600,AspectRatio->1];

Show[Table[Graphics[sampP],{2000}],Graphics[{PointSize->.01,Point[{x0[[1]],x0[[2]]}]}],ImageSize->600]

Cond[x]
Show[Table[Graphics[sampP],{2000}],Graphics[{PointSize->.01,Point[{x0[[1]],x0[[2]]}]}],gr,ImageSize->600]

Uszkodzenie mózgu

Wizualizacja graficzna uzyskanej sieci

p0=Table[{0,1.i/ni-1/2./ni},{i,1,ni}];
p1=Table[{1,1.i/n1-1/2./n1},{i,1,n1}];
p2=Table[{2,1.i/n2-1/2./n2},{i,1,n2}];
p3=Table[{3,1.i/no-1/2./no},{i,1,no}];

lines1=Table[Graphics[{If[W1[[i,j]]>0,Red,Blue],Thickness[Abs[W1[[i,j]]]/800],Line[{p0[[i]],p1[[j]]}]}],{i,1,ni},{j,1,n1}];

lines2=Table[Graphics[{If[W2[[i,j]]>0,Red,Blue],Thickness[Abs[W2[[i,j]]]/800],Line[{p1[[i]],p2[[j]]}]}],{i,1,n1},{j,1,n2}];

lines3=Table[Graphics[{If[W3[[i,j]]>0,Red,Blue],Thickness[Abs[W3[[i,j]]]/800],Line[{p2[[i]],p3[[j]]}]}],{i,1,n2},{j,1,no}];


Show[Graphics[PointSize[.02]],lines1,lines2,lines3,Graphics[Point[p0]],Graphics[Point[p1]],Graphics[Point[p2]],Graphics[Point[p3]],ImageSize->900,AspectRatio->.4]

Teraz wróć wyżej i ...

Sieć dla problemu "koła w kwadracie" uzyskana wcześniej jest inna!

Podsumowanie

Uwagi podsumowujące działanie wielowarstwowego perceptronu:

- Długie uczenie na dużej próbce polepsza wyniki.
- Użycie tej samej próbki możana powtarzać w iteracjach.
- Duże sieci, z wieloma warstwami i neuronami, mają wieksze możliwości.
- Wynikowe sieci “robiące tę samą czynność” mogą mieć zupełnie różne zestawy wag (lokalne minima), ale działają z podobną efektywnością. 
- Krawędzie o małych wagach można przeciąć, nie zaburzając zbytnio diałania sieci (uszkodzenie mózgu).
- Problem lokalnych minimów może być pokonany poprzez inne metody ustalania optymalnych wag.

Zauważmy też, że zaimplementowane powyżej rozpoznawanie figur nie ma nic wspólnego z geometrią i jest czysto “neuronowe”.

Dlaczego to działa? 

Pojedynczy neuron MCP z ostrym progiem: (Subscript[w, 0])+(Subscript[w, 1])(Subscript[x, 1])+..+(Subscript[w, n]) (Subscript[x, n]) = 0 jest równaniem hiperpłaszczyzny, zatem możemy oddzielić punkty należące do dwóch kategorii: nad i pod hiperpłaszczyzną. W ogólności wieloneuronowy perceptron składa takie warunki, będąc w stanie rozpoznać bardziej skomplikowane kategorie (obszary).

Show[Import["/home/broniows/local/Computer_Methods/neuron/ambr.png","PNG"],ImageSize->700]
(* Krystyna Ambroch, "Sztuczne sieci neuronowe", Gdańsk *)

ZADANIE 7

0) Przeczytaj notatki “Neural Networks” (Christos Stergiou, Dimitrios Siganos),  plik nn.pdf.
1) Sprawdź, czy sieć neuronowa z uczeniem poprzez propagację wsteczną podana na wykładzie może rozpoznać inne figury, a także figury z dziurami w środku, np. pierścień w kwadracie, itd. Może zaistnieć konieczność dodania kolejnych warstw do sieci.
2) [dla fizyków] Przelicz wyprowadznie algorytmu propagacji wstecznej na podstawie notatek back_prop.pdf.

Podsumowanie

Uwagi podsumowujące działanie wielowarstwowego perceptronu:

- Długie uczenie na dużej próbce polepsza wyniki.
- Użycie tej samej próbki możana powtarzać w iteracjach.
- Duże sieci, z wieloma warstwami i neuronami, mają wieksze możliwości.
- Wynikowe sieci “robiące tę samą czynność” mogą mieć zupełnie różne zestawy wag (lokalne minima), ale działają z podobną efektywnością. 
- Krawędzie o małych wagach można przeciąć, nie zaburzając zbytnio diałania sieci (uszkodzenie mózgu).
- Problem lokalnych minimów może być pokonany poprzez inne metody ustalania optymalnych wag.

Zauważmy też, że zaimplementowane powyżej rozpoznawanie figur nie ma nic wspólnego z geometrią i jest czysto “neuronowe”.

Dlaczego to działa? 

Pojedynczy neuron MCP z ostrym progiem: (Subscript[w, 0])+(Subscript[w, 1])(Subscript[x, 1])+..+(Subscript[w, n]) (Subscript[x, n]) = 0 jest równaniem hiperpłaszczyzny, zatem możemy oddzielić punkty należące do dwóch kategorii: nad i pod hiperpłaszczyzną. W ogólności wieloneuronowy perceptron składa takie warunki, będąc w stanie rozpoznać bardziej skomplikowane kategorie (obszary).

Show[Import["/home/broniows/local/Computer_Methods/neuron/ambr.png","PNG"],ImageSize->700]
(* Krystyna Ambroch, "Sztuczne sieci neuronowe", Gdańsk *)

Wykład 8 
Sieci neuronowe III

Oko kraba podkowiastego / tłumienie poboczne
Sieci Hopfielda

Pamiętacie złudzenie Macha?

(* funkcja intensywności *)
low=0.2; up=0.8;xl=300; xu=600;
p[x_]:=low /; x< xl
p[x_]:=up/; x>xu
p[x_]:= low+(x-xl)/(xu-xl)*(up-low) /;x ≥xl && x≤ xu
pic=Transpose[Table[{p[i]},{i,0,900}]];
ArrayPlot[pic,PlotRange->{0,1},ImageSize->700,AspectRatio->1,Frame->False]

Plot[p[x],{x,0,1000}]

Nasz brat krab podkowiasty pozwoli je nam zrozumieć!

Show[Import["/home/broniows/kielce/Computer_Methods_2015/neuron/Limule.jpg","JPG"],ImageSize->500]

Haldan Keffer Hartline (December 22, 1903 – March 17, 1983) was an American physiologist who was a cowinner (with George Wald and Ragnar Granit) of the 1967 Nobel Prize in Physiology or Medicine for his work in analyzing the neurophysiological mechanisms of vision. Hartline investigated the electrical responses of the retinas of certain arthropods, vertebrates, and mollusks because their visual systems are much simpler than those of humans and are thus easier to study. He concentrated his studies on the eye of the horseshoe crab (Limulus polyphemus). Using minute electrodes in his experiments, he obtained the first record of the electrical impulses sent by a single optic nerve fibre when the receptors connected to it are stimulated by light. He found that the photoreceptor cells in the eye are interconnected in such a way that when one is stimulated, others nearby are depressed, thus enhancing the contrast in light patterns and sharpening the perception of shapes. Hartline thus built up a detailed understanding of the workings of individual photoreceptors and nerve fibres in the retina, and he showed how simple retinal mechanisms constitute vital steps in the integration of visual information [Wikipedia].

Sieć Hopfielda

Show[Import["/home/broniows/kielce/Computer_Methods_2015/neuron/Hopfield-net.png","PNG"],ImageSize->300]

Model ze sprzężeniem zwrotnym dla N neuronów z tej samej warstwy:

(Subscript[f, n+1])=(Subscript[f, n])+ϵ (s+W.(Subscript[f, n]) - (Subscript[f, n]))

(Subscript[f, i]) - sygnał po iteracji i (wektor o wymiarze równym liczbie neuronów w warstwie, N)
s - sygnał (wektor o wymiarze N)
ϵ - parametr
W - symetryczna macierz sprzężenia o wymiarach N x N

Dla stanu ustalonego (po wielu iteracjach) mamy (Subscript[f, n]) ⟶ f, zatem
 
 s+W.f - f=0,
 f=(((1-W))^(-1)) s.         (***)

Przyjmiemy przykładowo 

(Subscript[W, ii])=0,
(Subscript[W, ij])=-a exp(-|i-j|/b) dla i≠j, a>0

tj. tłumienie jest najsilniejsze dla bliskich sasiadów i spada z odległością z charakterystyczna skalą b.

ns=30;       (* liczba neuronów N *)
b=6;           (* parametr spadku tłumienia z odległoscią neuronów i, j *)
a=0.2;     (* siła tłumienia *)
iter=100;  (* liczba iteracji *)
eps=.1;    (* parametr "szybkości uczenia " *)

(* przykładowy sygnał - paradoks Macha *)
s=Join[Table[0,{i,ns/3}],Table[3. i/ns,{i,ns/3}],
Table[1,{i,ns/3}]];

g0=ListPlot[s,PlotRange->{{0,30},{-0.1,1.1}},ImageSize->400,AxesLabel->{StyleForm["x",FontSize->18],StyleForm["sygnał wejściowy",FontSize->18]}]

f=Table[.001Random[],{ns}]; (* poczatkowy stan *)
W=Table[-a Exp[-Abs[i-j]/b],{i,ns},{j,ns}]; (* wagi *)

(* neuron nie tłumi sam siebie *)
Do[W[[i,i]]=0,{i,1,ns}]

ListDensityPlot[W,InterpolationOrder->0,Frame->True,FrameLabel->{StyleForm["i",FontSize->18],StyleForm["j",FontSize->18]},PlotLabel->StyleForm["Wagi sprzężenia",FontSize->18]]

ListPlot[Inverse[IdentityMatrix[ns]-W].s,ImageSize->400,AxesLabel->{StyleForm["x",FontSize->18],StyleForm["stan ustalony",FontSize->18]}]

(* implementacja jednej iteracji *)
T[f_]:=f+eps (s+W.f-f);

git=ListPlot[Table[Nest[T,f,i],{i,1,iter}],Joined->True,PlotRange->{{0,30},{-.2,1.0}},ImageSize->600,AxesLabel->{StyleForm["x",FontSize->18],StyleForm["odpowiedź w kolejnych iteracjach",FontSize->18]}]

rescale=2.5; (* aby lepiej porównac rysunki *)
g1=ListPlot[rescale Nest[T,f,iter],Joined->True,PlotRange->{{0,30},{-.4,1.3}}];

Show[g0,g1,Graphics[Text[StyleForm[iter "iterations",FontSize->18],{ns/2,-0.4}]],PlotRange->All,ImageSize->600,AxesLabel->{StyleForm["x",FontSize->18],StyleForm["sygnał wejściowy i przetworzony",FontSize->18]}]

Widać, że w modelu oka kraba płaskie fragmenty sygnału wejściowego są odwzorowywane ze zniekształceniem. Jakościowo odpowiada to złudzeniu Macha. W szczególności, dla x=10 mamy minimum (najjaśniejszy “pasek”), a dla x=20 (lokalne) maksimum (ciemny pasek). Tak więc prosty model ujemnego sprzężenia zwrotnego między receptorami siatkówki wystarcza do jakościowego wyjaśnienia zjawiska.

Po co Natura zadaje sobie tyle trudu, a przy okazji zmysły nas “oszukują? Jest tak, ponieważ ten sam mechanizm przetwarzania obrazu prowadzi do zwiększenia kontrastu/ostrości widzenia.

(* inny przykładowy sygnał *)
s=Table[3^2/((x-15.)^2+3^2),{x,1,30}];

g0=ListPlot[s,PlotRange->{{0,30},{-0.1,1.1}},ImageSize->400,AxesLabel->{StyleForm["x",FontSize->18],StyleForm["sygnał wejściowy",FontSize->18]}]

f=Table[.001 Random[],{ns}];

rescale=1.6;
g1=ListPlot[rescale Nest[T,f,iter],Joined->True,PlotRange->{{0,30},{0,1.0}}];

Show[g0,g1,Graphics[Text[StyleForm[iter"iterations",FontSize->18],{ns/2,-0.4}]],PlotRange->All,ImageSize->600,AxesLabel->{StyleForm["x",FontSize->18],StyleForm["sygnał wejściowy i przetworzony",FontSize->18]}]

Uzyskany obraz (linia ciagła) jest węższy od oryginału (kropki), sieć dokonała więc użytecznej korekcji/wyostrzenia obrazu.

Komentarze

Dla sieci Hopfielda (z symetrycznymi wagami (Subscript[w, ji])=(Subscript[w, ij ])) można wprowadzić pojęcie energii

E=-((1 )/2)(∑_ij) (Subscript[f, i]) (Subscript[w, ij]) (Subscript[f, j]) +(∑_i)(Subscript[f, i]) (Subscript[s, i]).

Warunki na minimum dają 

((∂E)/(∂f))=-(∑_j) (Subscript[w, ij]) (Subscript[f, j]) +(Subscript[s, i])=0

(ponadto macierz w musi być dodatnio określona), skąd f=(w^(-1))s (jest to ten sam warunek, co (***) powyżej dla w=1-W). A zatem sieci Hopfielda są równoważne modelowi spinowemu ze sprzeżeniami (Subscript[w, ij]) w zewnętrznym polu s! Wielkość (Subscript[s, i]) można też interpretować jako próg neuronu MCP.

Ponieważ powyższy warunek ma postać macierzową f=M s, podpada pod dyskutowane wcześniej modele pamieci heteroasocjatywnej i autoasocjatywnej (różnica tkwi jedynie w odwracaniu macierzy i związanym z tym “dostępie” poprzez iteracje. Tym samym sieci Hopfielda służą do modelowania pamięci neuronowych i w tym sensie są równoważne perceptronowi.

Sieć (maszyna) Boltzmanna

Jest to uogólnienie symetrycznej sieci Hopfielda dla (Subscript[f, i])=0,1 (dwa stany), wprowadzające termalny szum statystyczny. Formalnie, statystyczna funkcja rozdziału

Z=(∑_a) (e^(-(Subscript[E, a])/kT)), 

gdzie a jest konfiguracją całego układu a (Subscript[E, a]) jego zdefiniowaną powyżej energią - patrz wykład o sieciach spinowych (Metropolis, wyżarzanie, ...). Wprowadza się też czasem węzły (neurony) widoczne i ukryte. Uczenie maszyn Boltzmanna napotyka na problemy dla większej liczby neuronów i silnych połaczeń. Związane jest to z termalną aktywacja neuronów. Pozytywnym aspektem jest fakt, że temperatura pozwala na wyjście z lokalnych minimów (wyżarzanie).

Uczenie nienadzorowane / programowanie dynamiczne

Omawiane we wcześniejszych wykładach uczenie nadzorowane (patrz perceptron) potrzebuje nauczyciela. Biologicznie nie jest to realistyczne, bo wielu czynności uczymy się “na bieżąco” bez zewnetrznego nadzoru. Formalnie problem uczenia nienadzorowanego wiąże się z klasyfikacją danych (podziałem na klastry), przy czym nie wiemy z góry, jakie te dane są i jak są podzielone.

Quit[]

Wyobraźmy sobie, że kosmita wylądował na Ziemi i uczy się rozpoznawać rasy w miarę, jak spotyka kolejnych przedstawicieli naszego gatunku. Kolor skóry możemy reprezentować jako trzywymiarowy wektor barw w standardzie RGB, tj. postaci (r,g,b), gdzie wszystkie współrzędne są z przedziału [0,1]. Ponieważ opisany poniżej algorytm wymaga wektorów znormalizowanych (o tej samej długości), wprowadzamy fikcyjny czwarty wymiar w nastepujacy sposób:

x=(r,g,b,±(Sqrt[3-(r^2)-(g^2)-(b^2)])),  x.x=3.

Poniżej tworzeni są przedstawiciele poszczególnych ras. Ich kolory skóry fluktuują, ale w dopuszczalnych dla danej rasy granicach.

murzyn:=Module[{p},p={Random[Real,{0,.1}],Random[Real,{0,.1}],Random[Real,{0,.1}]};{p[[1]],p[[2]],p[[3]],Sqrt[3-p[[1]]^2-p[[2]]^2-p[[3]]^2]}]

bialy:=Module[{p},p={Random[Real,{.9,1}],Random[Real,{.9,1}],Random[Real,{.9,1}]};{p[[1]],p[[2]],p[[3]],-Sqrt[3-p[[1]]^2-p[[2]]^2-p[[3]]^2]}]

indianin:=Module[{p},p={Random[Real,{.9,1}],Random[Real,{.5,.5}],Random[Real,{.4,.5}]};{p[[1]],p[[2]],p[[3]],Sqrt[3-p[[1]]^2-p[[2]]^2-p[[3]]^2]}]

chinczyk:=Module[{p},p={Random[Real,{.9,1}],Random[Real,{.9,1}],Random[Real,{.45,.5}]};{p[[1]],p[[2]],p[[3]],Sqrt[3-p[[1]]^2-p[[2]]^2-p[[3]]^2]}]

(* próbka, której bedziemy dalej używać *)
tlumek=RandomSample[Join[Table[murzyn,{18}],Table[bialy,{20}],Table[indianin,{13}],Table[chinczyk,{22}]]];

(* próbka rezerwowa, która działa! *)
tlumek={{0.06296005388553759`,0.02516664198136456`,0.07983047325916161`,1.7288810737830784`},{0.06849824620533464`,0.06440964830488469`,0.024665785124328432`,1.7293209611047438`},{0.982161564111472`,0.9157861194184073`,0.9045965396563762`,0.6151418908710345`},{0.9316135518447106`,0.9598053361014183`,0.9836611993438572`,0.49323478356248424`},{0.9861731076152291`,0.5`,0.4183047096601366`,1.2658924803049707`},{0.9700434542074078`,0.5`,0.4642980386605903`,1.2623165325088987`},{0.9601775798855734`,0.5`,0.45279686747999864`,1.2739835210415336`},{0.9792262579544475`,0.992180542139847`,0.4784575957600452`,-0.9098197824795449`},{0.9898368135944611`,0.954059615973083`,0.41235006888068537`,-0.9695157308246997`}};

(* i-ta osoba z próbki z odciętą czwartą współrzędną *)
person[i_]:=Module[{},pp=Delete[tlumek[[i]],4];{RGBColor[pp[[1]],pp[[2]],pp[[3]]],PointSize[Large],Point[pp]}]

(* wizualizacja próbki w przestrzeni RGB *)
pl=Show[Graphics3D[Table[person[i],{i,1,Length[tlumek]}]]]

Teraz tworzymy (rekursywną) sieć jednowarstwową o imax=4 neuronach i inicjalizujemy losowo  wektory wagowe (znormalizowane do 1). Neuronów powinno być tyle, ile kategorii (klastrów danych).

imax=4;

(* inicjalizacja losowa *)
Do[w[i]={1-.2Random[],Random[],1-.5Random[],1-2.Random[]};no[i]=Sqrt[w[i].w[i]];w[i]=w[i]/no[i],{i,1,imax}]

(* inicjalizacja z pomocą osobników *)
w[1]=tlumek[[1]]/Sqrt[3];
w[2]=tlumek[[2]]/Sqrt[3];
w[3]=tlumek[[3]]/Sqrt[3];
w[4]=tlumek[[4]]/Sqrt[3];

Show[pl,Table[Graphics3D[{Opacity[0.25],Sphere[Delete[w[i],4] Sqrt[3],.04]}],{i,1,imax}]]

Algorytm jest nastepujący: bierzemy osobnika s, dla każdego znajdujemy wektor wagowy (Subscript[w, i]), który ma z nim najwiekszy iloczyn skalarny. Nastepnie obracamy (Subscript[w, i]) w kierunku s (możliwe są tu różne przepisy). Powtarzamy procedure dla wszystkich osobników i ew. powtarzamy całość nauki kilkukrotnie. W wyniku procedury  otrzymujemy wektory reprezentatywne dla każdej z klas. 

Przyjęta poniżej strategia nazywa się “
winner takes all”, tj. obracamy tylko jeden wektor najbliższy danemu osobnikowi. Uogólnieniem jest strategia “winner takes most”, która umożliwia uaktualnianie  również dalszych wektorów wagowych.

eps=.1;

Do[s=tlumek[[j]]/Sqrt[3];
tt=Table[w[i].s,{i,1,imax}];
iup=Position[tt,Max[tt]][[1]][[1]];
wup=w[iup]+eps(s-w[iup]);
nup=Sqrt[wup.wup];
w[iup]=wup/nup,{j,1,Length[tlumek]}]

Show[pl,Table[Graphics3D[{Opacity[0.2],Sphere[Delete[w[i],4] Sqrt[3],.04]}],{i,1,imax}],Table[Graphics3D[{Text[i,Delete[w[i],4] Sqrt[3]]}],{i,1,imax}]]

Nowi napotkani osobnicy klasyfikowani są wg istniejacych wektorów reprezentatywnych - są najbliżej wektora, który reprezentuje ich rasę.

proxy[a_,b_]:=Sqrt[Sum[(a[[i]]/Sqrt[3]-b[[i]])^2,{i,1,3}]]

mu=murzyn;
bi=bialy;
in=indianin;
ch=chinczyk;

{proxy[mu,w[1]],proxy[mu,w[2]],proxy[mu,w[3]],proxy[mu,w[4]]}
{proxy[bi,w[1]],proxy[bi,w[2]],proxy[bi,w[3]],proxy[bi,w[4]]}
{proxy[in,w[1]],proxy[in,w[2]],proxy[in,w[3]],proxy[in,w[4]]}
{proxy[ch,w[1]],proxy[ch,w[2]],proxy[ch,w[3]],proxy[ch,w[4]]}

Sieci Kohonena

[rozdział oparty o książkę Raúl Rojas, Neural Networks, a Systematic Introduction (neuron_book.pdf)]

Show[Import["/home/broniows/local/Computer_Methods_2014/neuron/map_vis.png","PNG"],ImageSize->900]

(Subscript[□, □])Samoorganizujące się odwzorowanie (fiński neurobiolog Teuvo Kohonen, 1970’) zachowujace topologię. 

 - plakietkowanie: jeśli input należy do plakietki, jeden neuron odpala
 - stwarzane jest “rusztowanie” przestrzeni, którą obrazujemy (b. ważne, podstawowe zadanie niezbędne do postrzegania)
 - sieci Kohonena tworzą mapy przestrzeni wejściowej w samoorganizujący się sposób (uczenie nienadzorowane)
 - liniowe lub planarne struktury neuronowe mózgu przetwarzają sygnały wielowymiarowe (percepcja jest wielowymiarowa - głębia, kolor, tekstura,...)
 - zachowanie topologii - obiekty bliskie (z tego samego otoczenia) pozostają po odwzorowaniu blisko
 - mapowanie w mózgu nie jest, lub nie jest do końca, zdeterminowane genetycznie

Show[Import["/home/broniows/local/Computer_Methods/neuron/koh1.png","PNG"],ImageSize->600]

Algorytm Kohonena

n-wymiarowe wagi (Subscript[w, 1]) , (Subscript[w, 2]) , . . . , (Subscript[w, m]) dla m neuronów są losowo inicjalizowane. Stała uczenia η oraz funcja sprzężenia sąsiadów φ, zależna od promienia sąsiedztwa, są odpowiednio dobrane.

1 : Wybierz (losowo) wektor inputowy x z przestrzeni inputowej z odpowiednim rozkładem prawdopodobieństwa.
2 : Wybieramy neuron k, dla którego euklidesowa odległość (Subscript[w, i ])od x jest najmniejsza.
3 : Wagi są uaktualniane wg przepisu

(Subscript[w, i]) ⟶ (Subscript[w, i]) + ηφ(i, k)(x − (Subscript[w, i])),  i = 1, . . . , m.

4 : Stop, jeśli osiągnieto maksymalną liczbę kroków, w przeciwnym razie zmodyfikuj η i φ zgodnie z planem i idź do 1.

(Uaktualnianie sąsiadów jest kluczowe dla zachowania topologii - bliskie neurony z dużym prawdopodobieństwem też odpalają)

Show[Import["/home/broniows/local/Computer_Methods_2014/neuron/koh3.png","PNG"],ImageSize->1300]

Show[Import["/home/broniows/local/Computer_Methods_2014/neuron/koh2.png","PNG"],ImageSize->900]

Powyżej wynik odwzorowania kwadratu w dwuwymiarową sieć Kohonena po 100, 1000, 5000 i 10000 iteracjach.

Podsumowanie sieci neuronowych

- model obliczania całkowicie inny od maszyny von Neumanna (automat komórkowy, komputer), masywnie równoległe przetwarzanie we wzajemnie ze sobą oddziałujących neuronach
- ta sama klasa obliczalności wg Churcha-Turinga
- neuron MCP (McCullogh-Pitts), wagi jako nośnik informacji (pamięć, program, dane, ...)
- perceptron (Rosenblatt) jako sieć jednokierunkowa, sieci rekursywne Hopfielda, Boltzmanna, Kohonena, ... - ustalanie wag wg różnych strategii
- model pamięci (heteroasocjatywna - skojarzenia, autoasocjatywna - naprawa utraconych/niekompletnych danych)
- uczenie nadzorowane (propagacja wsteczna) i nienadzorowane (dynamiczne)
- algorytm klasyfikacji (klasteryzacji) danych w uczeniu nienadzorowanym - systemy diagnostyczne, detekcja trudnych do wychwycenia korelacji
- problemy związane ze zbieżnością, lokalnymi minimami - cała dziedzina neuromatematyki
- analogia do modeli spinowych (fizyka statystyczna, wyżarzanie)
- droga do zrozumienia działania mózgu (świadomość, inteligencja, kreatywność, wolna wola)
- jednoczesny niesłychany rozwój techniki obrazowania mózgu (Human Connectome Project, ...) - możliwość coraz lepszego porównania modeli z realizacją biologiczną
- sztuczna inteligencja ...

Show[Import["/home/broniows/local/Computer_Methods_2014/neuron/con.jpg","JPG"],ImageSize->700]

ZADANIE 8

0) Przeczytaj pierwszą część rozdziału 15 z książki Raúl Rojas, Neural Networks, a Systematic Introduction (neuron_book.pdf).

1) Wytłumacz jakościowo złudzenie Hermanna w oparciu o tłumienie poboczne. Co można wywnioskować o działaniu receptorów plamki żółtej i pozostałych receptorów siatkówki?

Do wyboru 2) lub 3):
2) Zaprogramuj odwzorowywanie kwadratu w dwuwymiarową sieć Kohonena (dającą efekt, jak na powyższym rysunku).
3) Zaprogramuj odwzorowywanie kwadratu w jednowymiarową sieć Kohonena.

Wykład 9
Dynamika nieliniowa i chaos I

m.in. na podstawie 
- Steven H. Strogatz, Nonlinear Dynamics and Chaos, Westview Press 1994 
- Predrag Cvitanović i in., Chaos Book, http://chaosbook.dk/

Quit[]

Równania różnicznowe nieliniowe rzędu 1
(tzw. teoria jakościowa rr)

Prosty przykład

Scałkujmy równanie 

dx/dt = sin(x), x((Subscript[t, 0]))=(Subscript[x, 0])
dt = dx/sin(x)

Quit[]

s[x_]=Integrate[1/Sin[x],x]

t[x_,x0_]=s[x]-s[x0]+t0

Zebranie logarytmów do kupy

rule={Log[a__]+Log[b__]->Log[a b],Log[a__]-Log[b__]->Log[a /b]}

so=t[x,x0]//.rule

Plot[Table[t[x,x0]/.t0->0,{x0,Pi/24,2Pi,Pi/12}],{x,-13,10},AxesLabel->{"x","t"}]

Okresowość x -> x + 2 π jest oczywista, wystarczy więc rozważać jedną gałąź

so

inv=Solve[so==t,x]

inv=Simplify[inv/.C[1]->k, {Im[t-t0]==0,Element[k,Integers]}]

Plot[Table[x/.inv/.k->0 /.t0->0,{x0,Pi/24,2Pi,Pi/12}],{t,-8,8},AxesLabel->{"t-(Subscript[t, 0])","x"}]

To samo z pomocą wbudowanej procedury

sol=DSolve[{x'[t]==Sin[x[t]],x[t0]==x0},x[t],t]

Trajektorie o bliskich warunkach początkowych (po tej samej stronie punktu x=0) pozostają blisko dla wszystkich czasów

Punkty stacjonarne (punkty stałe)

Są to punkty, dla których "przepływ" dx/dt = 0. Dla Naszego przykładu warunkiem tym jest sin(x)=0. Prowadzi to do punktów stałych stabilnych i niestabilnych. Przepływ oznaczony strzałkami prowadzi od puntktów niestabilnych do stabilnych (czas zawsze biegnie do przodu).

Show[Plot[Sin[x],{x,-10,10},AspectRatio->.3],Graphics[Table[Disk[{Pi+2k Pi,0},{.25,.1}],{k,-2,1}]],Graphics[Table[Circle[{2k Pi,0},{.25,.1}],{k,-2,1}]],Graphics[Table[Arrow[{{1+2k Pi,0},{1.8+2k Pi,0}}],{k,-2,1}]],Graphics[Table[Arrow[{{1+2k Pi,0},{-1.8+2k Pi,0}}],{k,-2,1}]],AxesLabel->{"x","dx/dt=sin(x)"}]

Simplify[Normal[Series[x'[t]==Sin[x[t]],{x[t],k Pi,1}]],Element[k,Integers]]

Równanie logistyczne (modelowanie populacji)
Verhulst 1938

(przeskalowanie zmiennej x i t - patrz poniżej)

eq=x'[t]==x[t](1-x[t])

Show[Plot[x(1-x),{x,0,1}],Graphics[Disk[{1,0},{.025,.01}]],Graphics[Circle[{0,0},{.025,.01}]],Graphics[Arrow[{{.3,0},{.5,0}}]],Graphics[Arrow[{{1.3,0},{1.2,0}}]],AxesLabel->{"x","dx/dt=x(1-x)"},PlotRange->All,AspectRatio->.55]

sol=DSolve[{eq,x[0]==x0},x[t],t][[1]]

pl=Plot[{x[t]/.sol/.x0->.2,x[t]/.sol/.x0->.6,x[t]/.sol/.x0->1.2,x[t]/.sol/.x0->-.001},{t,0,6},PlotStyle->{{Red,Thickness[.01]}},PlotRange->{-.3,1.5},AxesLabel->{"t","x"}]

Przepływ (dx, dt) w przestrzeni (x, t)

plv=VectorDensityPlot[{1,x(1-x)},{t,0,6},{x,-.3,1.3},FrameLabel->{"t","x"}]

Show[plv,pl]

Dla problemu populacji ujemne rozwiązania są formalne

Linearyzacja wokół punktów stałych - analiza lokalnej stabilności

Rozwijamy równanie wokół punktu stałego (Subscript[x, 0]). Niech x=(Subscript[x, 0])+δ, wówczas

dx/dt = dδ/dt = f((Subscript[x, 0])+δ) = f((Subscript[x, 0]))+f'((Subscript[x, 0])) δ  + ... = λ δ + ...

skąd δ(t) = c (e^(λt)). Tak więc znak pochodnej (znak λ) determinuje stabilność, a wielkość określa “prędkość” zmierzania do lub oddalania się od punktu stałego. Możliwa jest sytuacja λ=0 (np. dx/dt=(x^2)). Odpowiada to punktowi siodłowemu. Stałe λ nazywamy wykładnikami Lyapunova, a powyższe warunki stanowią kryterium stabilności Lyapunova.

Uwagi

Rozważamy układy n równań autonomicznych, tj. takich, gdzie funkcja f nie zależy od t. Równania nieautonomiczne można sprowadzić do autonomicznych dodając dodatkową zmienną (Subscript[x, n+1])=t i równanie (Subscript[dx, n+1])/dt = 1.

Tw. Jeśli f i f’ są ciągłe na przedziale A zawierajacym (Subscript[x, 0]), to istnieje przedział B zawierający (Subscript[x, 0]) na którym istnieje jednoznaczne rozwiązanie problemu Cauchy’ego dx/dt = f,  x((Subscript[t, 0]))=(Subscript[x, 0]). Kontrprzykładem jest f(x)=(x^(1/3)) (rozwiązania osobliwe, brak jednoznaczności). 

Jednoznaczność oznacza, że różne krzywe całkowe nie mają punktów wspólnych.

Niemożliwy jest ruch drgający dla równań autonomicznych rzędu pierwszego na R (przepływ jest w kierunku punktów stabilnych).

Quit[]

ss=DSolve[{x'[t]==x[t]^(1/3),x[0]==0},x[t],t]//FullSimplify

Plot[{0,x[t]/.ss},{t,0,3},Frame->True,FrameLabel->{"t","x"}]

Widać, że dla t=0 mamy dwa rozwiązania o x(0)=0.

Bifurkacje (rozdwojenia)

Zajmiemy się teraz zależnością rozwiązań nieliniowych rr rzędu pierwszego od zewnętrznych parametrów (teoria katastrof, Rene Thom).

Bifukacja punktu siodłowego

Rozważmy równanie dx/dt = a+(x^2)

GraphicsGrid[{{Show[Plot[-1+x^2,{x,-1.2,1.2},PlotRange->{-1,2},PlotLabel->"a=-1"],Graphics[Disk[{-1,0},{.08,.16}]],Graphics[Circle[{1,0},{.08,.16}]]],Show[Plot[0+x^2,{x,-1.2,1.2},PlotRange->{-1,2},PlotLabel->"a=0"],Graphics[{Gray,Disk[{0,0},{.09,.18}]}]],Show[Plot[1+x^2,{x,-1.2,1.2},PlotRange->{-1,2},PlotLabel->"a=1"]]}},ImageSize->800]

Plot[{-Sqrt[-a],Sqrt[-a]},{a,-3,0},PlotStyle->{Blue,Dashed},PlotRange->{{-3,1},Automatic},Frame->True,AxesOrigin->{1,-1.6},FrameLabel->{StyleForm["a",FontSize->18],StyleForm["punkty stałe",FontSize->18]},PlotLabel->StyleForm["bifurkacja punktu siodłowego",FontSize->18]]

g[x_,a_]=x[t]/.DSolve[{x'[t]==a+x[t]^2},x[t],t][[1]]

Zmiana charakteru rozwiązania w punkcie bifurkacji a=0.

g[x_,1] =x[t]/.DSolve[{x'[t] ==  1+x[t]^2},x[t],t][[1]]
g[x_,0] =x[t]/.DSolve[{x'[t] ==  0+x[t]^2},x[t],t][[1]]
g[x_,-1]=x[t]/.DSolve[{x'[t]==-1+x[t]^2},x[t],t][[1]]

Bifukacja transkrytyczna

Rozważmy równanie dx/dt = a x-(x^2)

GraphicsGrid[{{Show[Plot[-x-x^2,{x,-1.2,1.2},PlotRange->{-2,1},PlotLabel->"a=-1"],Graphics[Disk[{-1,0},{.08,.16}]],Graphics[Circle[{0,0},{.08,.16}]]],Show[Plot[-x^2,{x,-1.2,1.2},PlotRange->{-2,1},PlotLabel->"a=0"],Graphics[{Gray,Disk[{0,0},{.09,.18}]}]],Show[Plot[x-x^2,{x,-1.2,1.2},PlotRange->{-2,1},PlotLabel->"a=1"],Graphics[Disk[{1,0},{.08,.16}]],Graphics[Circle[{0,0},{.08,.16}]]]}},ImageSize->800]

Show[Plot[{0,a},{a,-2,0},PlotStyle->{Blue,Dashed},Frame->True,FrameLabel->{StyleForm["a",FontSize->18],StyleForm["punkty stałe",FontSize->18]},PlotLabel->StyleForm["bifurkacja transkrytyczna",FontSize->18]],Plot[{a,0},{a,0,2},PlotStyle->{Blue,Dashed},Frame->True,FrameLabel->{StyleForm["a",FontSize->18],StyleForm["punkty stałe",FontSize->18]}],PlotRange->All,AxesOrigin->{-2,-2}]

h[x_,a_]=x[t]/.DSolve[{x'[t]==a x[t]-x[t]^2},x[t],t][[1]]

Zmiana charakteru rozwiązania

h[x,1]
Limit[h[x,a],a->0]
h[x,-1]

Bifukacja widelcowa

Rozważmy równanie dx/dt = a x-(x^3)

GraphicsGrid[{{Show[Plot[-x-x^3,{x,-1.2,1.2},PlotRange->{-2,2},PlotLabel->"a=-1"],Graphics[Disk[{0,0},{.08,.16}]]],Show[Plot[-x^3,{x,-1.2,1.2},PlotRange->{-2,2},PlotLabel->"a=0"],Graphics[{Gray,Disk[{0,0},{.09,.18}]}]],Show[Plot[x-x^3,{x,-1.2,1.2},PlotRange->{-2,2},PlotLabel->"a=1"],Graphics[Disk[{1,0},{.08,.16}]],Graphics[Disk[{-1,0},{.08,.16}]],Graphics[Circle[{0,0},{.08,.16}]]]}},ImageSize->800]

Show[Plot[{0},{a,-1.5,0},PlotStyle->{Blue},Frame->True,FrameLabel->{StyleForm["a",FontSize->18],StyleForm["punkty stałe",FontSize->18]},PlotLabel->StyleForm["bifurkacja widelcowa nadkrytyczna",FontSize->18]],Plot[{Sqrt[a],-Sqrt[a],0},{a,0,2},PlotStyle->{Blue,Blue,Dashed},Frame->True,FrameLabel->{StyleForm["a",FontSize->18],StyleForm["punkty stałe",FontSize->18]}],PlotRange->{-2,2},AxesOrigin->{-1.5,-2}]

“Nadkrytyczna” oznacza, że przy zmianie parametru punkt stabilny rozszczepia się na 2 punty stabilne i jeden niestabilny. 

Występuje też bifurkacja podkrytyczna, o równaniu prototypowym x'=ax+(x^3), gdzie punkt stabilny “anihiluje” z dwoma niestabilnymi, tworząc punkt niestabilny.W wyniku tego punkt stabilny znika i system musi radykalnie zmienić zachowanie, skacząc do innego punktu stabilnego lub do nieskończoności. Poniżej obrazek.

Show[Plot[{0},{a,0,1.5},PlotStyle->{Dashed},Frame->True,FrameLabel->{StyleForm["a",FontSize->18],StyleForm["punkty stałe",FontSize->18]},PlotLabel->StyleForm["bifurkacja widelcowa podkrytyczna",FontSize->18]],Plot[{Sqrt[-a],-Sqrt[-a],0},{a,-2,0},PlotStyle->{{Dashed,Blue},{Dashed,Blue},Blue},Frame->True,FrameLabel->{StyleForm["a",FontSize->18],StyleForm["punkty stałe",FontSize->18]}],PlotRange->{-2,2},AxesOrigin->{-2,-2}]

Ilustracja dla bifurkacji widelcowej nadkrytycznej:

k[x_,a_]=x[t]/.DSolve[{x'[t]== a x[t]-x[t]^3},x[t],t][[1]]

Zmiana charakteru rozwiązania w punkcie bifurkacji:

k[x,1]=x[t]/.DSolve[{x'[t]== x[t]-x[t]^3},x[t],t][[1]]
k[x,0]=x[t]/.DSolve[{x'[t]== -x[t]^3},x[t],t][[1]]
k[x,-1]=x[t]/.DSolve[{x'[t]== -x[t]-x[t]^3},x[t],t][[1]]

Przykład: inwazja zwójkowatych, czyli katastrofa wierzchołka

n(t) - populacja owadów
R - stała określająca tempo rozmnażania się owadów
K - stała określająca populację nasycenia owadów przy braku ptaków
p(t) - tempo, w jakim ubywa owadów z powodu żerowania ptaków
A, B - stałe określające żerowanie

eq=n'[t]==R n[t](1-n[t]/K)-p[t]
p[t]=B n[t]^2/(A^2 +n[t]^2)

Poprzez przekslowanie t i n(t) można pozbyć sie 2 z 4 stałych.

eq1=eq /.{n[t]->A x[t],n'[t]->B  x'[t],t->τ}

eq1[[1]]/B == eq1[[2]]/B  /.K-> k A /.R-> B/A r//FullSimplify

eq2=x'[t]==r x[t](1-x[t]/k) - x[t]^2/(1+x[t]^2)

Powyższe równanie zależy od dwóch stałych, bedących kombinacjami początkowych stałych R, K, A i B:

k=K/A, r=AR/B.

Trywialnym punktem stałym jest x(t)=0. Analiza graficzna pozwala zrozumieć ulokowanie pozostałych punktów stałych.

aa=eq2[[2]]/. x[t]->y

lewa=aa[[1]]/y
prawa=-aa[[2]]/y

Przecięcia ukazują punkty stabilne. Jest ich, w zalezności od wartości r i k od 0 do 3.

Show[Plot[{(lewa/.{r->.45,k->9.5}),prawa},{y,0,10},AxesLabel->{StyleForm["y",FontSize->18],""}],Graphics[Text[StyleForm["k",FontSize->16],{9.7,.04}]],Graphics[Text[StyleForm["r",FontSize->16],{0.4,.47}]]]

Bifurkacja zachodzi wówczas, gdy niebieska prosta jest styczna do żółtej krzywej. Wtedy ich pochodne w punkcie styczności są sobie równe.

(* równość pochodnych *)
eqs=D[lewa,y]==D[prawa,y]//Simplify

so=Solve[eqs,r]

(* równość funkcji *)
solk=Solve[lewa==prawa/.so[[1]],k][[1]]

Zauważmy, że k > 0  implikuje  y > 1

solr=Solve[lewa==prawa/.solk,r]

Tak więc bifurkacja punktu siodłowego zachodzi w przestrzeni (k,r) wzdłuż krzywej parametrycznej o równaniu:

poi[y_]={k,r}/.solk/.solr

ParametricPlot[poi[y],{y,1.04,12},AspectRatio->1,PlotRange->{{0,25},{0,.7}},AxesLabel->{StyleForm["k",FontSize->18],StyleForm["r",FontSize->18]}]

Rozwiązanie "brutalne"

y3d=y/.Solve[lewa==prawa //Simplify,y]

Plot3D[{y3d[[1]], y3d[[2]], y3d[[3]]}, {k, 2, 20}, {r, 0.1, 0.7}, Mesh -> None, 
  PlotPoints -> 35, AxesLabel -> {StyleForm["k", FontSize -> 18], 
    StyleForm["r", FontSize -> 18], 
    StyleForm["y", FontSize -> 18]}, 
  PlotLabel -> StyleForm["punkty stale populacji zwójkowatych",FontSize->16]]

(żółty i niebieski - punkty stałe stabilne, zielony - punkty stałe niestabilne)

Teoria katastrof - klasyfikacja zachowań układów nieliniowych

Teoria katastrof zajmuje się systemami dynamicznymi, w szczególności sytuacjami, gdy ciągła zmiana parametru kontrolnego prowadzi do jakościowyh zmian zachowania systemu. Jej autorem jest francuski matematyk René Thom, lata 50. XX w. Klasyfikację można rozszerzyć dla większej liczby parametrów i na równania wyższego rzędu.

Vladimir Arnold:

    A0 - a non-singular point: V = x.
    A1 - a local extremum, either a stable minimum or unstable maximum 
                                                                                            V = +/- (x^2) + a x.
    A2 - the fold
    A3 - the cusp
    A4 - the swallowtail
    A5 - the butterfly
    Ak - a representative of an infinite sequence of one variable forms V=x^{k+1}+ ...
    D4- - the elliptical umbilic
    D4+ - the hyperbolic umbilic
    D5 - the parabolic umbilic
    Dk - a representative of an infinite sequence of further umbilic forms
    E6 - the symbolic umbilic V = x^3+y^4+a x y^2 +bxy+cx+dy+ey^2
    E7
    E8

Zadanie 9

0) Zapoznaj się z książką “Nonlinear dynamics and chaos”, Steven H. Strogatz i przeczytaj interesujący Cię fragment dot. teorii katastrof w zastosowaniu do któregoś z omawianych tam modeli.

1) Zanalizuj w podobny sposób jak na wykładzie model populacji ryb z członem “odławiania”:

n’(t)=R n(t)(1-n(t)/K) - A n(t)/(B+n(t)), R,K,A,B>0.

2) Rozważ następujący model epidemii (Kermack, McKendrick, 1927):

x(t) - populacja ludzi zdrowych
y(t) - populacja ludzi chorych
z(t) - zmarli

x’=-κ x y,
y’= κ x y - λ y,
z’= λ y.

a) Przedyskutuj sens powyższych równań.
b) Użyj pierwszego i ostatniego równania aby pokazać, że x’/x=-κ/λ z’, a po scałkowaniu x(t)=(Subscript[x, 0]) exp[-κ/λ z(t)].
c) Pokaż, że x+y+z=n=const.
d) W oparciu o b) i c) pokaż, że 

z'=λ(n-z-(Subscript[x, 0])  exp[-κ/λ z(t)]).

e) Zanalizuj to równanie przy pomocy naszych metod.

Wykład 10
Dynamika nieliniowa i chaos II

Rząd 2

linia → płaszczyzna

Punkt stały jest na płaszczyźnie
Analiza stabilności: linearyzacja wokół punktu stałego, 
x=u (e^(λt)),  Au = λ u → równanie własne dla macierzy dwuwymiarowej. 
Jeśli dla wszystkich λ zachodzi Re(λ)<0, to x jest stabilny. 
Możliwość orbit zamkniętych dla przestrzeni (R^2) - periodyczność.
Cykle graniczne (zjawisko stricte nieliniowe).

Oscylator harmoniczny

Quit[]

ho={v'[t]==-w^2 x[t],x'[t]==v[t],x[0]==x0,v[0]==0}

s={x[t],v[t]}/.DSolve[ho,{x[t],v[t]},t][[1]]/.w->1

Show[VectorDensityPlot[{v,-x},{x,-1,1},{v,-1,1}],ParametricPlot[Table[s,{x0,0,1,.2}],{t,0,2Pi},PlotStyle->{Red,Thickness[.01]}],FrameLabel->{"x","v"}]

(środek)

A={{0,-w^2},{1,0}}

Eigenvalues[A]

Eigensystem[A]//MatrixForm

Środek w (0,0)

Układy zachowawcze, energia

ContourPlot[1/2 (x^2+v^2),{x,-1,1},{v,-1,1},Contours->{0.01^2,.1^2,.2^2,.3^2,.4^2,.5^2,.6^2,.7^2}]

Oscylator anharmoniczny

pl=ContourPlot[1/2 (x^4-x^2+v^2),{x,-1.3,1.3},{v,-1,1},Contours->{-.34^2,-.3^2,-.2^2,-.1^2,0,.1^2,.2^2,.3^2,.4^2,.5^2,.6^2,.7^2},FrameLabel->{StyleForm["x",FontSize->18],StyleForm["v",FontSize->18]}]

Dwa środki i siodło

Show[VectorDensityPlot[{v,-2x^3+x},{x,-.9,.9},{v,-.5,.5}],FrameLabel->{"x","v"}]

Ogólna klasyfikacja dla Ax = λx

(* S. H. Strogatz *)
Show[Import["/home/broniows/kielce/Computer_Methods_2015/chaos/2p.png","PNG"],ImageSize->900]

Ogólna sytuacja
Możliwe cykle graniczne - zjawisko wymagające nieliniowości!

(* S. H. Strogatz *)
Import["/home/broniows/kielce/Computer_Methods_2015/chaos/p3.png","PNG",ImageSize->700]

Trywialny przykład dający cykl graniczny r=1 we współrzędnych biegunowych:

dr/dt=r(1-r)
dϕ/dt=1

Oscylator van der Pola

Przykład nietrywialnego cyklu granicznego, który istnieje dla μ > 0
(sieci z lampami elektronowymi).

Quit[]

eq=x''[t]+μ(x[t]^2-1)x'[t]+x[t]==0

tmax=40;

vdp=NDSolve[{eq/.μ->1.7,x[0]==.15,x'[0]==.5},x[t],{t,0,tmax}][[1]];

pl=Plot[x[t]/.vdp,{t,0,tmax},AxesLabel->{StyleForm["t",FontSize->18],StyleForm["x(t)",FontSize->18]}]

nxp[t_]=D[x[t]/.vdp,t];

pl2=ParametricPlot[{x[t]/.vdp,nxp[t]},{t,0,tmax},AxesLabel->{StyleForm["x(t)",FontSize->18],StyleForm["v(t)",FontSize->18]}]

GraphicsGrid[{{pl,pl2}},ImageSize->700]

Definicje stabilności

1) ∃δ>0:|(SuperStar[x])-x(0)|<δ ⟹ (lim_(t → ∞))x(t)=(SuperStar[x])  - punkt przyciągajacy

2) ∀ϵ>0 ∃δ>0:|(SuperStar[x])-x(0)|<δ ∧ t>0  ⟹ |x(t)-(SuperStar[x])|<ϵ  - stabilny w sensie Lyapunova

1) i 2) = asymptotycznie stabilny


Np. środek jest 2) ale nie 1).

Tw. Poincarego - Bendixsona

Jeśli trajektoria T układu 2-wymiarowego o f(x) klasy (Subscript[C, 1]) zawarta jest dla  t>(Subscript[t, 0 ])w ograniczonym obszarze, a obszar ten nie zawiera punktów stałych, to T dąży do cyklu granicznego.

Bifurkacja Hopfa

{{b,2},{-1,0.8}}//MatrixForm

Eigenvalues[{{b,2},{-1,0.8}}]

Manipulate[t=Eigenvalues[{{b,2},{-1,0.8}}];ListPlot[{{Re[t[[1]]],Im[t[[1]]]},{Re[t[[2]]],Im[t[[2]]]}},PlotStyle->PointSize[Large],PlotRange->{{-2,1},{1.5,-1.5}}],{b,-2.65,1}]

Część urojona wartości własnych powoduje oscylacje. Po destabilizacji punktu staje się on niestabilną spiralą. Jeśli ruch jest ograniczony, to układ dąży na mocy twierdzenia P-B do cyklu granicznego. 

Taką sytuację nazywamy nadkrytyczną bifurkacją Hopfa: stabilna spirala przechodzi w spiralę niestabilną otoczoną stabilnym cyklem granicznym. Promień cyklu zachowuje się jak (((a - (Subscript[a, crit])))^(1/2)). 

Przykład:

r'=a r-(r^3),
ϕ'=ω+b (r^2).

Manipulate[Plot[a r-r^3,{r,0,1},AxesLabel->{StyleForm["r",FontSize->18],StyleForm["r'",FontSize->18]},PlotRange->{-.3,.3}],{a,-.4,2}]

Podkrytyczna bifurkacja Hopfa zachodzi w sytuacji, gdy punkt stabilny otoczony jest granicznym cyklem niestabilnym, który z kolei jest otoczony stabilnym cyklem granicznym. Dla krytycznej wartości parametru cykl niestabilny "anihiluje" punkt stabilny. Układ musi wiec dążyć do (odległego) stabilnego cyklu granicznego.

Przykład:

r'=a r+(r^3)-(r^5),
ϕ'=ω+b (r^2).

Manipulate[Plot[a r+r^3-r^5,{r,0,1},AxesLabel->{StyleForm["r",FontSize->18],StyleForm["r'",FontSize->18]},PlotRange->{-.3,.3}],{a,-.4,2}]

Import["/home/broniows/kielce/Computer_Methods_2015/chaos/p3.png","PNG",ImageSize->700]

Niestabilna spirala z B otoczona cyklem granicznym D.

Więcej wymiarów i więcej nieliniowości...

Show[Import["/home/broniows/kielce/Computer_Methods_2015/chaos/tab.png","PNG"],ImageSize->1200]

Mapa Poincarego

Import["/home/broniows/kielce/Computer_Methods_2015/chaos/pomap.png","PNG",ImageSize->500]

(Subscript[x, n+1])=P((Subscript[x, n]))

Cykle zmienione w punkty, prostszy opis.

Równania Lorenza

Edward Lorenz (1963) - bardzo uproszczony model cieplnej konwekcji atmosferycznej

Quit[]

Import["/home/broniows/kielce/Computer_Methods_2015/chaos/lor.png","PNG",ImageSize->500]

Grzanie oceaniczne od dołu, chłodzenie od góry, pole grawitacyjne, lepkość, przewodnictwo cieplne, 3 wymiary zredukowane do dwóch (symetria translacyjna), wymuszony ruch w “rurze” jak na rysunku.

x ~ q  - bezwymiarowa prędkość
y ~ (Subscript[T, 3]) - bezwymiarowa różnica temperatur między prądami wznoszącymi i opadającymi 
z ~ (Subscript[T, 2])-(Subscript[T, 1  ])- bezwymiarowe odstępstwo temperatur od równowagi
t - bezwymiarowy czas
   
    σ - liczba Prandtla - charakteryzująca lepkość ośrodka
    ρ - liczba Rayleigha - charakteryzująca przewodnictwo cieplne
    β - stała charakteryzująca rozmiary obszaru
    
    x'==σ (y-x),
y'==x (ρ-z)-y,
z'==x y-β z.

σ=10;
ρ=28;
β=8/3;

lor={x'[t]==σ (y[t]-x[t]),y'[t]==x[t](ρ-z[t])-y[t],z'[t]==x[t]y[t]-β z[t]}

tmax=65;
z0=23;
ls=NDSolve[{lor,x[0]==1,y[0]==1,z[0]==z0},{x[t],y[t],z[t]},{t,0,tmax}][[1]];

pl=ParametricPlot3D[{x[t],y[t],z[t]}/.ls,{t,0,tmax},PlotRange->All,ColorFunction->Function[{x,y,z, t}, Hue[2t]],PlotPoints->200,ImageSize->600,AxesLabel->{StyleForm["x(t)",FontSize->18],StyleForm["y(t)",FontSize->18],StyleForm["z(t)",FontSize->18]}]

Nieregularne (aperiodyczne/chaotyczne) oscylacje, trajektoria ograniczona.

Plot[x[t]/.ls,{t,0,tmax},PlotRange->All,PlotStyle->Blue,PlotPoints->200,ImageSize->800,AxesLabel->{StyleForm["t",FontSize->18],StyleForm["x(t)",FontSize->18]},AspectRatio->.2]

Efekt motyla

eps=0.00001
ls2=NDSolve[{lor,x[0]==1,y[0]==1,z[0]==z0+eps},{x[t],y[t],z[t]},{t,0,tmax}][[1]];

ParametricPlot3D[({x[t],y[t],z[t]}/.ls)-({x[t],y[t],z[t]}/.ls2),{t,0,tmax},PlotRange->All,PlotStyle->Blue,PlotPoints->200,AxesLabel->{StyleForm["Δx(t)",FontSize->18],StyleForm["Δy(t)",FontSize->18],StyleForm["Δz(t)",FontSize->18]},ImageSize->600]

LogPlot[Sqrt[Apply[Plus,(({x[t],y[t],z[t]}/.ls)-({x[t],y[t],z[t]}/.ls2))^2]],{t,0,tmax},PlotRange->All,AxesLabel->{StyleForm["t",FontSize->18],StyleForm["Δ(t)",FontSize->18]},ImageSize->800,AspectRatio->.2]

Po pewnym czasie pojawia się “chaotyczna” odległość dwóch poczatkowo bardzo bliskich trajektorii.

Definicja chaosu

Chaos = aperiodyczne zachowanie układu deterministycznego (ograniczonego) w długim czasie, przejawiające czułość na warunki początkowe.

Mieszanie topologiczne - po pewnym czasie każdy obszar przestrzeni fazowej będzie "zmieszany" z dowolnym innym obszarem (zob. "wałkowanie ciasta" poniżej).

Gęste orbity periodyczne - przez dowolnie małe otoczenie każdego punktu przestrzeni fazowej przechodzi orbita periodyczna (niestabilna).

Struktura fraktalna, samopodobieństwo przy zmianie skali.

Import["/home/broniows/kielce/Computer_Methods_2015/chaos/ciasto.png","PNG",ImageSize->600]

Dziwny atraktor

Atraktor (punkt stabilny, cykl graniczny stabilny)

Podzbiór A przestrzeni fazowej spełniający warunki:
 1) x((Subscript[t, 0])) ϵ A ⇒ x(t) ϵ A 
 2) Istnieje zbiór otwarty U zawierający A taki, że jeśli x((Subscript[t, 0])) ϵ U, to w granicy t → ∞ odległość x(t) od A dąży do zera.
 3) A jest minimalny, tj. nie zawiera podzbioru właściwego spełniającego 1) i 2).

Dziwny atraktor - atraktor, który przejawia czułość na warunki początkowe, tzn. dwa punkty będące blisko siebie na dziwnym atraktorze w czasie (Subscript[t, 0  ])oddalają się od siebie eksponencjalnie z upływem czasu. 

Wszystie znane dziwne atraktory mają strukturę fraktalną.

Mapa Lorenza

fu[t_]=z[t]/.ls;

Plot[fu[t],{t,0,tmax},PlotRange->All,PlotStyle->Blue,PlotPoints->200,ImageSize->800,AxesLabel->{StyleForm["t",FontSize->18],StyleForm["z(t)",FontSize->18]},AspectRatio->.2]

Znajdziemy wszystkie maksima z (t) pomiędzy (Subscript[t, 0]) a (Subscript[t, 1])

allz[t0_,t1_,st_]:=Module[{t00,t11,tmin,tmax,ru},
ze={};t00=t0;t11=t00+st/8;tmin=t00-st/2;tmax=t00+st/2;
While[t11<t1,ru=FindMaximum[fu[y],{y,t00,t11,tmin,tmax}, MaxIterations->200];
zz=y/.ru[[2]];
If[tmin<zz && zz<=tmax,ze=Append[ze,ru[[1]]]];
t00=t00+st;t11=t11+st;tmin=tmin+st;tmax=tmax+st]]

allz[0,65,.02]

ze

tze=Table[{ze[[i]],ze[[i+1]]},{i,1,Length[ze]-1}]

Show[ListPlot[tze],Plot[x,{x,32,42}],AxesLabel->{StyleForm["(Subscript[x, n])",FontSize->18],StyleForm["(Subscript[x, n + 1])",FontSize->18]},PlotLabel->StyleForm["mapa Lorenza",FontSize->18],ImageSize->450]

Pewna prostota wyłania się z chaosu!

Lorenz: Ponieważ |f’(x)|>1, każda orbita zamknięta jest niestabilna (dowód poprzez rozwinięcie wokół punktu stałego mapy (f^p)(x) - Strogatz).

Zadanie 10

0) Wyprowadzenie równań Lorenza znajduje się w pliku Lorenz_derivation.pdf
1) Rozważ (z pomocą Mathematiki) harmonicznie wymuszony tłumiony oscylator nieliniowy,

x''(t)=-(x^3)+ x + γ x’(t) + A sin(ωt)

dla różnych wartości parametrów i wykaż, że układ może być chaotyczny.
2) Skonstruuj dla powyższego układu mapę Poincarego, zdefiniowaną jako zbiór punktów (x(t),y(t)) dla ωt = 2πn, gdzie n jest całkowite.

Wykład 11
Dynamika nieliniowa i chaos III

Analiza przebiegów czasowych

Quit[]

xs[t_]=x[t]/.NDSolve[{x''[t]==x[t]-x[t]^3+ Sin[3 t],x[0]==0.1,x'[0]==1},x[t],{t,0,100}][[1]]

Plot[xs[t-5]-xs[t],{t,5,100}]

Proste kryterium: brak okresowości x(t+δ)-x(t) w przypadku układów chaotycznych

Mapa logistyczna

Mapy są często prostsze w analizie od równań różniczkowych

(Subscript[x, n+1])=f((Subscript[x, n]))

(Poincare, Lorenz)

Quit[]

f[x_,r_]= r x(1-x)

rr=.8

rr=2.9

rr=3.2

rr=3.5

rr=3.65

rr

pl=Plot[{x,f[x,rr]},{x,0,1}];
x0=Random[]
x1=f[x0,rr]
pl=Show[pl,Graphics[Line[{{x0,0},{x0,x1},{x1,x1}}]]]

step:=Module[{},x0=x1;
x1=f[x0,rr];
pl=Show[pl,Graphics[Line[{{x0,x0},{x0,x1},{x1,x1}}]]]]

step

pl=Plot[{x,f[x,rr]},{x,0,1}];

tt=Table[step,{300}];

ListAnimate[tt,AnimationRepetitions->1]

Tw. Banacha o punkcie stałym

Def. Funkcja f: X→X spełnia warunek Lipschitza (jest zwężeniem), jeśli d(f(x),f(y)) ≤ q d(x,y), gdzie 0<q<1 jest stałą Lipschitza.

Tw. (Banach) Zwężenie na przestrzeni metrycznej zupełnej posiada jeden i tylko jeden punkt stały (SuperStar[x])= f((SuperStar[x])).

“Mapa Polski rzucona na ziemię ma dokładnie jeden punkt, który sie pokryje z odpowiadającym mu punktem na ziemi.” 

Uwaga: Cykle o okresie p można rozważać jako punkty stałe odwzorowania (f^n)(x).

Clear[x]
rr=2.9;
xc=x/.Solve[x==f[x,rr],x][[2]];

q=.99999;
Show[ContourPlot[Abs[f[x,rr]-f[y,rr]]/Abs[x-y]-q,{x,.9 xc-.01,1.101xc+.01},{y,.9 xc-.01,1.1xc+.01},Contours->{0}],Graphics[{Red,Disk[{xc,xc},.005 xc+.0005]}]]

Mapa logistyczna nie jest zwężeniem w całej dziedzinie

Bifurkacje zdwajania okresu dla mapy logistycznej

Clear[x,r]

x1=x/.Solve[x==f[x,r],x]

(* nowe zera *)
eq=(x-f[f[x,r],r])/x/(x-x1[[2]])

eq=eq//FullSimplify

x2=x/.Solve[eq==0,x]

x2/.r->3.01

rr=3.3
Plot[{x,f[f[x,rr],rr]},{x,0,1}]

rr=3.49
Plot[{x,f[f[f[f[x,rr],rr],rr],rr]},{x,0,1}]

Clear[x,r]

x2

eq=(x-f[f[f[f[x,r],r],r],r])/x/(x-x1[[2]])  //FullSimplify

ss=Solve[eq==0,x]

ss/. r->2.+Sqrt[6]

Numeryczne badanie cykli mapy logistycznej

Quit[]

f[x_,r_]= r x(1-x)

tt={};

Timing[Do[
x=Random[];
Do[x=f[x,r],{200}];
Do[tt=AppendTo[tt,{r,x=f[x,r]}],{100}],{r,2.8,4,.0021}]]

ListPlot[tt,PlotStyle->PointSize[.001],ImageSize->900,AxesLabel->{StyleForm["r",FontSize->18],StyleForm["(SuperStar[x])",FontSize->18]}]

Widzimy podwajania okresu

Samopodobieństwo

tt={};

Timing[Do[
x=Random[];
Do[x=f[x,r],{400}];
Do[tt=AppendTo[tt,{r,x=f[x,r]}],{100}],{r,3.5,3.6,.0002}]]

ListPlot[tt,PlotStyle->PointSize[.001],ImageSize->600,PlotRange->{.32,.43},AxesLabel->{StyleForm["r",FontSize->18],StyleForm["(SuperStar[x])",FontSize->18]}]

Samopodobieństwo!

Intermitencja

Obszar przejścia miedzy regularnością a chaosem

r0=3.8282
r0=4
x=.8;

tt=Table[{i,x=f[x,r0]},{i,1,200}];
ListPlot[{tt,tt},PlotStyle->PointSize[.005],AspectRatio->.2,ImageSize->800,Joined->{True,False},AxesLabel->{StyleForm["n",FontSize->18],StyleForm["(Subscript[x, n])",FontSize->18]}]

tt=Table[{i,x=f[x,r0]},{i,201,400}];
ListPlot[{tt,tt},PlotStyle->PointSize[.005],AspectRatio->.2,ImageSize->800,Joined->{True,False},AxesLabel->{StyleForm["n",FontSize->18],StyleForm["(Subscript[x, n])",FontSize->18]}]

tt=Table[{i,x=f[x,r0]},{i,401,600}];
ListPlot[{tt,tt},PlotStyle->PointSize[.005],AspectRatio->.2,ImageSize->800,Joined->{True,False},AxesLabel->{StyleForm["n",FontSize->18],StyleForm["(Subscript[x, n])",FontSize->18]}]

Stałe Feygenbauma [Mitchell Feygenbaum (1944-)]

Moduł znajdujący punkt nieciągłości funkcji

Quit[]

newt[f_,x0_,x1_,eps_]:=Module[{xa,xb},xa=x0;xb=x1; While[(xa-xb)^2>eps^2  && f[xa]≠f[xb],xc=(xa+xb)/2;If[f[xa]==f[xc],xa=xc,xb=xc]];xc]

Sprawdzenie dla funkcji schodkowej Heavyside' a

newt[UnitStep,-1,.3,.00001]

Znajdowanie pierwiastków n - krotnego złożenia mapy logistycznej

Manipulate[
g[x_]=r x(1-x);
x/.NSolve[x==Nest[g,x,n],x],{n,1,4,1},{r,1,4.1,.001}]

gr zlicza rzeczywste dodatnie pierwiastki funkcji h(r, x) w zmiennej x dla ustalonego r

gr[r_]:=Module[{},
tt=x/.NSolve[x==h[r,x],x,Reals];rr=0;
Do[If[tt[[i]]>0,rr++],{i,1,Length[tt]}];rr]

Badamy teraz (2^n)-krotnie iterowane mapy logistyczne

Mapa iterowana raz

h[r_,x_]=Nest[r # (1-#)&,x,1]
Plot[{x,h[.9,x]},{x,0,1}]
Plot[{x,h[1.3,x]},{x,0,1}]
Plot[gr[r],{r,0,2}]
Print["n=1, pierwiastek pojawia się dla r=",newt[gr,.9,1.1,.000001]]

h[r_,x_]=Nest[r # (1-#)&,x,2]
Plot[{x,h[2.9,x]},{x,0,1}]
Plot[{x,h[3.1,x]},{x,0,1}]
Plot[gr[r],{r,2.8,3.2}]
Print["n=2, pierwiastek pojawia się dla r=",newt[gr,2.8,3.18,.000001]]

h[r_,x_]=Nest[r # (1-#)&,x,4]
Plot[{x,h[3.44,x]},{x,0,1}]
Plot[{x,h[3.47,x]},{x,0,1}]
Plot[gr[r],{r,3.44,3.46}]
Print["n=4=(2^2), pierwiastek pojawia się dla r=",newt[gr,3.44,3.46,.00000001]]

h[r_,x_]=Nest[r # (1-#)&,x,8];
Plot[{x,h[3.54,x]},{x,0,1}]
Plot[{x,h[3.55,x]},{x,0,1}]

(* liczy się zbyt długo! *)
(* Plot[gr[r],{r,3.541,3.55}]
newt[gr,3.541,3.55,.00000001] *)

Kolejne wartości r, dla których okres mapy logistycznej się podwaja :

fg={3.0, 3.44949,3.54409,3.5644,3.568759}
fglim=3.569946

Stosunki Feigenbauma :

Do[Print[(fg[[i+1]]-fg[[i]])/(fg[[i+2]]-fg[[i+1]])],{i,1,3}]

Limit:  δ = 4.669201609... - stała “Przyrody” dot. chaosu

Układ dochodzi do chaosu w sposób uniwersalny!

(* S. H. Strogatz *)
Show[Import["/home/broniows/kielce/Computer_Methods_2015/chaos/figa.png","PNG"],ImageSize->500]

α=(Subscript[d, n])/(Subscript[d, n+1]) → -2.502907875095892822283902873218... (pierwsza stała)
δ=(Subscript[Δ, n])/(Subscript[Δ, n+1])→ 4.669201609... (druga stała Feigenbauma)

Uniwersalność Feigenbauma zachodzi dla wszystkich map unimodalnych, tj. funkcji różniczkowalnych wypukłych ku górze z pojedynczym kwadratowym maksimum

Pomiary fizyczne stałej δ

Show[Import["/home/broniows/local/Computer_Methods/chaos/f_ex.png","PNG"],ImageSize->800]

Fraktale

Zbiór Julii

Gaston Julia (1893–1978)

Quit[]

f[z_]=z^2-1.01+.411I
out=100
test[z_]:=Abs[z]<3
JuliaCount[z0_]:=Length[NestWhileList[f,N[z0],test,1,out]];

JuliaCount[0]
JuliaCount[1+.5I]

DensityPlot[JuliaCount[x+I y],{x,-1.9,1.9},{y,-1.2,1.2},PlotPoints->300,AspectRatio->1.2/1.9,ImageSize->600]

Dla funkcji zespolonych całkowitych zbiór Julii jest granicą między obszarem z którego mapa f "wyprowadza" (Subscript[z, n]) do nieskończoności i obszarem, w którym (Subscript[z, n ])pozostaje ograniczone.

Zbiór Mandelbrota

Benoit Mandelbrot (1924, Warszawa - 2010)

Teraz pytamy, dla jakich wartości liczby zespolonej C mapa (Subscript[x, n+1])=((Subscript[x, n])^2)+C, (Subscript[x, 0])=0, pozostawia orbitę ograniczoną. Taki zbiór liczb zespolonych C to właśnie zbiór Mandelbrota.

Quit[]

clowerleft =2.0 +I 2.0;
cupperright =2.0 +I 2.0;
cstep =0.05;
rmax =2.0;
Nmax =50;

mandelbrotEquation= Compile[{{C,_Complex},{Nmax,_Integer},{rmax,_Real}},Length[FixedPointList[#^2+ C&,C,Nmax,SameTest ->(Abs[#2]>  rmax&)]]];

mandelbrotEquation[0.5,100,3]

DensityPlot[mandelbrotEquation[x+y I,Nmax,rmax],{x,-2.3,1},{y,-1.2,1.2},PlotPoints->200,AspectRatio->.7,ImageSize->1200]

Stosunki promieni kolejnych "kół" dążą do δ Feigenbauma!
Samopodobieństwo, popkultura

Własności: fraktal, ograniczony (|c|<2), spójny
Pole: 1.50659177 ± 0.00000008
Wzdłuż osi rzeczywistej odpowiedniość z mapą logistyczną - idąc po osi rzeczywistej od prawej do lewej napotykamy na kolejne bifurkacje na granicach podobszarów, przy czym c=r/2 (1 - r/2)

Poszczególne podobszary (tzw. hiperboliczne) odpowiadają cyklom o długości n (n=1 - punkt stały)

Show[Import["/home/broniows/kielce/Computer_Methods_2015/chaos/mand_num.png","PNG"],ImageSize->450]

Zadanie 11

1) Pobaw się fraktalami Julii i Mandelbrota
2) Powtórz symulacje wykładu dla mapy eliptycznej 
	xn+1 →r  (Sqrt[1-((((Subscript[x, n])-1))^2)]), (Subscript[x, n]) ϵ [0,2], r ϵ [0,2],
	
3) oraz dla mapy trójkątnej
	xn+1 →r  (Subscript[x, n]) , (Subscript[x, n]) ϵ [0,1],
	xn+1 →r  (2-(Subscript[x, n])) , (Subscript[x, n]) ϵ [1,2],
	r ϵ [1,2].
	
4) Pokaż, że dla pewnej rzeczywistej stałej c i przetransformowanej zmiennej (Subscript[x, n ])mapę Mandelbrota można utożsamić z mapą logistyczną.

Wykład 12
Chaos kwantowy

Równanie Schroedingera jest liniowe → nie ma chaosu w ewolucji funkcji falowej w sensie wykładniczej czułości na warunki początkowe. Obowiązuje jednak zasada odpowiedniości: dla dużych liczb kwantowych (lub formalnie gdy stała Plancka dąży do zera) ewolucja kwantowa dąży do ewolucji klasycznej, w związku z tym pewne cechy klasycznego układu powinny być widoczne dla układu kwantowego o bardzo dużych liczbach kwantowych.

Chaos kwantowy = jakie własności ma układ kwantowy, którego klasyczny hamiltonian prowadzi do chaosu?

Przykład układu regularnego - bilard w prostokącie

Quit[]

Wymiary prostokąta są niewspółmierne

a=1
b=Pi/3

Równanie liniii przechodzącej przez ((Subscript[x, 0]),(Subscript[y, 0])) o współczynniku kierunkowym ((Subscript[vx, 0]),(Subscript[vy, 0]))

line[{x_,y_},{x0_,y0_},{vx0_,vy0_}]=-vy0 x+vx0 y+x0 vy0-y0 vx0

line[{x,y},{0,1},{1,0}]

Przecięcie linii z obrzeżem prostokąta

pl[{x0_,y0_},{vx0_,vy0_}]:=NSolve[{line[{x,y},{x0,y0},{vx0,vy0}]==0,x==a || x==0 || y==b || y==0,x>=0,x≤a,y>=0,y≤b},{x,y}]

{x,y}/.pl[{0,1},{1,0}]

Pozbywamy się z rozwiązania punktu wyjściowego

(* poi[{x0_,y0_},{vx0_,vy0_}]:=DeleteCases[{x,y}/.pl[{x0,y0},{vx0,vy0}],{x0,y0}][[1]] *)

poi[{x0_,y0_},{vx0_,vy0_}]:=Module[{kk,nn},kk={x,y}/.pl[{x0,y0},{vx0,vy0}];If[(kk[[1]][[1]]-x0)^2+(kk[[1]][[2]]-y0)^2<.001^2,nn=2];If[(kk[[2]][[1]]-x0)^2+(kk[[2]][[2]]-y0)^2<.001^2,nn=1];kk[[nn]]]

poi[{0,1},{1,0}]

Prędkość po odbiciu od ściany

vn[{vx0_,vy0_},{x1_,y1_}]:= Which[y1==0 || y1==b,{-vx0,vy0},x1==0 || x1==a,{vx0,-vy0}]

Nowy punkt odbicia i predkość

new[{{x0_,y0_},{vx0_,vy0_}}]:= {poi[{x0,y0},{vx0,vy0}],vn[{vx0,vy0},poi[{x0,y0},{vx0,vy0}]]}

Tworzymy teraz listę punktów odbicia

new[{{0.,.29501},{1.,.97}}]

n=new[{{0.,.2501},{1.,.711}}]
tt={n[[1]]}

Do[n=new[n];AppendTo[tt,n[[1]]],{700}]

Show[Graphics[{Blue,Thickness[.01],Line[{{0,0},{0,b},{a,b},{a,0},{0,0}}]}],Graphics[Line[tt]],ImageSize->900]

Ten sam problem potraktowany kwantowo

Quit[]

Widmo dla nieskończonej prostokątnej studni potencjału

a=1;
b=Pi/3;
(* b=1 *)

e[n_,m_]=Pi^2 (n^2/a^2 +m^2/b^2)
max=350;
em=Pi^2 max^2/Max[a^2,b^2];

ft:=Module[{},tt={};Do[Do[If[e[n,m]<em,tt=AppendTo[tt,e[n,m]//N]],{n,1,max}],{m,1,max}];tt=Sort[tt];]

ft
td=Table[tt[[i+1]]-tt[[i]],{i,1,Length[tt]-1}];
av=Mean[td]
td=td/av;
Length[td]

Show[Histogram[td,200,"PDF"],Plot[Exp[-s],{s,0,7},PlotRange->All,PlotStyle->Thickness[.005]],PlotRange->{{0,5},All},ImageSize->500,Frame->True,FrameLabel->{StyleForm["s",FontSize->16],StyleForm["P(s)",FontSize->16]}]

Rozkład odległości między poziomami w (genercznym) układzie klasycznie regularnym jest poissonowski, P(s)=(e^(-s))
(hipoteza Berry’ego-Tabora, 1977)

Bilard stadionowy klasyczny

Leonid Bunimovich

Quit[]

a=1
b=Pi/3

Równanie liniii przechodzącej przez ((Subscript[x, 0]),(Subscript[y, 0])) o współczynniku kierunkowym ((Subscript[vx, 0]),(Subscript[vy, 0]))

line[{x_,y_},{x0_,y0_},{vx0_,vy0_}]=-vy0 x+vx0 y+x0 vy0-y0 vx0

Przecięcie linii z obrzeżem stadionu

pls1[{x0_,y0_},{vx0_,vy0_}]:=NSolve[{line[{x,y},{x0,y0},{vx0,vy0}]==0,x^2+(y-b/2)^2==b^2/4,x<0},{x,y}]
pls2[{x0_,y0_},{vx0_,vy0_}]:=NSolve[{line[{x,y},{x0,y0},{vx0,vy0}]==0,(x-a)^2+(y-b/2)^2==b^2/4,x>a},{x,y}]
pls3[{x0_,y0_},{vx0_,vy0_}]:=NSolve[{line[{x,y},{x0,y0},{vx0,vy0}]==0,y==0,x<=a,x>=0},{x,y}]
pls4[{x0_,y0_},{vx0_,vy0_}]:=NSolve[{line[{x,y},{x0,y0},{vx0,vy0}]==0,y==b,x<=a,x>=0},{x,y}]

vyy=2

pls1[{0,0},{1,vyy}]
pls2[{0,0},{1,vyy}]
pls3[{0,0},{1,vyy}]
pls4[{0,0},{1,vyy}]

Kolejny punkt odbicia i prędkość po odbiciu

news[{{x0_,y0_},{vx0_,vy0_}}]:=Module[{au,xa,ya,nx,ny,ss},
au=pls1[{x0,y0},{vx0,vy0}];If[au≠{},{xa,ya}={x,y}/.au[[1]];If[(xa-x0)^2+(ya-y0)^2>.001^2,{x1,y1}={xa,ya};{nx,ny}={x1,y1-b/2}/Sqrt[x1^2+(y1-b/2)^2];ss=vx0 nx+vy0 ny;{vx1,vy1}={vx0-2ss nx,vy0-2ss ny}]];
au=pls2[{x0,y0},{vx0,vy0}];If[au≠{},{xa,ya}={x,y}/.au[[1]];If[(xa-x0)^2+(ya-y0)^2>.001^2,{x1,y1}={xa,ya};{nx,ny}={x1-a,y1-b/2}/Sqrt[(x1-a)^2+(y1-b/2)^2];ss=vx0 nx+vy0 ny;{vx1,vy1}={vx0-2ss nx,vy0-2ss ny}]];
au=pls3[{x0,y0},{vx0,vy0}];If[au≠{},{xa,ya}={x,y}/.au[[1]];If[(xa-x0)^2+(ya-y0)^2>.001^2,{x1,y1}={xa,ya};{nx,ny}={0,1};ss=vx0 nx+vy0 ny;{vx1,vy1}={vx0-2ss nx,vy0-2ss ny}]];
au=pls4[{x0,y0},{vx0,vy0}];If[au≠{},{xa,ya}={x,y}/.au[[1]];If[(xa-x0)^2+(ya-y0)^2>.001^2,{x1,y1}={xa,ya};{nx,ny}={0,-1};ss=vx0 nx+vy0 ny;{vx1,vy1}={vx0-2ss nx,vy0-2ss ny}]];
au=pls1[{x0,y0},{vx0,vy0}];If[au≠{},{xa,ya}={x,y}/.au[[2]];If[(xa-x0)^2+(ya-y0)^2>.001^2,{x1,y1}={xa,ya};{nx,ny}={x1,y1-b/2}/Sqrt[x1^2+(y1-b/2)^2];ss=vx0 nx+vy0 ny;{vx1,vy1}={vx0-2ss nx,vy0-2ss ny}]];
au=pls2[{x0,y0},{vx0,vy0}];If[au≠{},{xa,ya}={x,y}/.au[[2]];If[(xa-x0)^2+(ya-y0)^2>.001^2,{x1,y1}={xa,ya};{nx,ny}={x1-a,y1-b/2}/Sqrt[(x1-a)^2+(y1-b/2)^2];ss=vx0 nx+vy0 ny;{vx1,vy1}={vx0-2ss nx,vy0-2ss ny}]];
 au=pls3[{x0,y0},{vx0,vy0}];If[au≠{},{xa,ya}={x,y}/.au[[2]];If[(xa-x0)^2+(ya-y0)^2>.001^2,{x1,y1}={xa,ya};{nx,ny}={0,1};ss=vx0 nx+vy0 ny;{vx1,vy1}={vx0-2ss nx,vy0-2ss ny}]];
au=pls4[{x0,y0},{vx0,vy0}];If[au≠{},{xa,ya}={x,y}/.au[[2]];If[(xa-x0)^2+(ya-y0)^2>.001^2,{x1,y1}={xa,ya};{nx,ny}={0,-1};ss=vx0 nx+vy0 ny;{vx1,vy1}={vx0-2ss nx,vy0-2ss ny}]]; 
{{x1,y1},{vx1,vy1}}]

Tworzymy teraz listę punktów odbicia

n={{a/2.,0},{.33,Random[]}}
tt={n[[1]]}

Do[n=news[n];AppendTo[tt,n[[1]]],{700}]

Show[Graphics[{Blue,Thickness[.005],Line[{{0,0},{a,0}}]}],Graphics[{Blue,Thickness[.005],Line[{{0,b},{a,b}}]}],Graphics[{Blue,Thickness[.005],Circle[{0,b/2},b/2,{Pi/2,3Pi/2}]}],Graphics[{Blue,Thickness[.005],Circle[{a,b/2},b/2,{-Pi/2,Pi/2}]}],Graphics[Line[tt]],ImageSize->700]

Układ jest chaotyczny

Przypadek kwantowy bilardu stadionowego

Kwadrat funkcji falowej dla bardzo dalekiego poziomu dla bilardu stadionowego:

Show[Import["/home/broniows/kielce/Computer_Methods_2015/chaos/stadium.png","PNG"],ImageSize->800]

Problem sprowadza się do rozwiązania równania Helmholza z warunkiem brzegowym na skomplikowanym konturze. Rysunki zaczerpnięto z literatury.

W przypadku chaotycznym (symetrycznym względem odbicia czasu i po wydzieleniu przez symetrie) rozkład odległości między poziomami jest bliski rozkładowi Wignera (patrz dalsza część wykładu)

P(s)=π s/2 (e^(-π (s^2)/4)).

Tak więc ważnym wyróżnikiem układów kwantowych, których Hamiltoniany prowadzą klasycznie do chaosu bądź regularności, jest rozkład odległości między poziomami. W przypadku chaotycznym poziomy się odpychają i rozkład jest bliski rozkładu Wignera. W przypadku regularnym (z generycznymi parametrami, np. niewspółmierne długości boków prostokąta w bilardzie) otrzymujemy rozkład Poissona. Jest to zadziwiająca (do dziś nie udowodniona) i uniwersalna cecha układów kwantowych!

Show[Import["/home/broniows/kielce/Computer_Methods_2015/chaos/stadwig.png","PNG"],ImageSize->400]

Doświadczalny (!) rozkład odległości dla stadionu

Miraż kwantowy

W okolicach klasycznych orbit periodycznych dla niektórych stanów widać bardzo wyraźne wzmocnienie funkcji falowej - spektakularne zjawisko mirażu kwantowego dla układów chaotycznych

Show[Import["/home/broniows/kielce/Computer_Methods_2015/chaos/mir.png","PNG"],ImageSize->700]

Atom wodoru w polu magnetycznym

Dla dostatecznie silnego pola magnetycznego klasycznie następuje przejście do chaosu

Show[Import["/home/broniows/kielce/Computer_Methods_2015/chaos/hyd1.png","PNG"],ImageSize->1000]
Show[Import["/home/broniows/kielce/Computer_Methods_2015/chaos/hyd2.png","PNG"],ImageSize->1000]

(slajdy zaczerpnięte z Dominique Delande, Quantum Chaos)

Macierze przypadkowe

J. Wishart (1928)

Fizyka jądrowa, jądro złożone - Eugene Wigner 

Przypuszczenie Wignera (1957) dot. rozkładu odległości między poziomami energetycznymi

Hipoteza Bohigasa, Giannoniego i Schmita (BGS, 1984): statystyka spektralna generycznego układu kwantowego, którego wersja klasyczna przejawia chaos, jest opisywalna w ramach teorii macierzy przypadkowych (wyjątek od uniwersalności: diamagnetyczny atom litu).

Generowanie liczb losowych z danym rozkładem :

RandomVariate[NormalDistribution[]]

Gauss orthogonal ensemble (GOE)

Układy z symetrią odbicia czasu

Quit[]

n=1000
A=Table[0,{n},{n}];
Do[A[[i,i]]=RandomVariate[NormalDistribution[0,1]],{i,1,n}]
Do[Do[A[[i,j]]=RandomVariate[NormalDistribution[0,1/Sqrt[2]]];A[[j,i]]=A[[i,j]],{i,1,j-1}],{j,1,n}]

A=A/Sqrt[2n];

Prawo półkola Wignera

tt=Sort[Eigenvalues[A]];
Show[Histogram[tt,n/40,"PDF"],Plot[2/Pi Sqrt[1-x^2],{x,-1,1}]]

Rozkład odległości poziomów jest bliski rozkładowi Wignera, ale nie identyczny. Poprawny (skomplikowany) rozkład podali Gaudin i Mehta (1960).

ss=Table[tt[[i+1]]-tt[[i]],{i,1,Length[tt]-1}];
ss=ss/Mean[ss];
Show[Histogram[ss,{.15},"PDF"],Plot[Pi s /2 Exp[-Pi s^2/4],{s,0,3}],PlotRange->{{0,3},All}]

Odpychanie poziomów! 
Bardzo dobrze widać jakościową różnicę z rozkładem Poissona, dyskutowanym powyżej dla przypadku regularnego (niechaotycznego).

Gauss unitary ensemble (GUE)

Układy bez symetrii odbicia czasu

n=1000
A=Table[0,{n},{n}];
Do[A[[i,i]]=RandomVariate[NormalDistribution[0,1/Sqrt[2]]],{i,1,n}]
Do[Do[A[[i,j]]=RandomVariate[NormalDistribution[0,1/2]]+I RandomVariate[NormalDistribution[0,1/2]];A[[j,i]]=Conjugate[A[[i,j]]],{i,1,j-1}],{j,1,n}]

A=A/Sqrt[2n];

tt=Sort[Eigenvalues[A]];
Show[Histogram[tt,n/20,"PDF"],Plot[2/Pi Sqrt[1-x^2],{x,-1,1}]]

Teraz rozkład odległości poziomów jest przybliżony funkcją 32(s^2)/(π^2) Exp(-4 (s^2)/π).

ss=Table[tt[[i+1]]-tt[[i]],{i,1,Length[tt]-1}];
ss=ss/Mean[ss];
Show[Histogram[ss,{0.15},"PDF"],Plot[32s^2/Pi^2 Exp[-4 s^2/Pi],{s,0,3}],PlotRange->{{0,3},All}]

Inne układy

Bilard Sinaia

Jeszcze jeden popularny przykład bilardu kwantowego: Bilard Sinaia, czyli “koło wycięte w kwadracie” (Yakov G. Sinai, 1963). Klasycznie układ jest chaotyczny i ergodyczny. Celem wydzielnia symetrii rozważa się układ w wycinku kwadratu, jak ukazano  w prawym górnym rogu poniższego rysunku.

Widmo doświadczalne (czerwony histogram) pochodzi z pomiaru wzbudzeń nadprzewodzącej wnęki rezonansowej o kształcie bilardu Sinaia.

Show[Import["/home/broniows/kielce/Computer_Methods_2015/chaos/sinai.png","PNG"],ImageSize->500]

Odległości między liczbami losowymi mają rozkład Poissona

Quit[]

Zabawa w zmianę “sąsiedniości” di

di=1
tt=Table[Random[],{50000}];
tt=Sort[tt];

ss=Table[tt[[i+di]]-tt[[i]],{i,1,Length[tt]-di}];
ss=ss/Mean[ss];
Show[Histogram[ss,{0.05},"PDF"],Plot[{Exp[-s],Pi s /2 Exp[-Pi s^2/4],32s^2/Pi^2 Exp[-4 s^2/Pi]},{s,0,4}],PlotRange->{{0,4},All}]

Różnica drugiego rzędu

ss2=Sort[Table[Abs[tt[[i+2di]]-2tt[[i+di]]+tt[[i]]],{i,1,Length[tt]-2di}]];
ss2=ss2/Mean[ss2];
Show[Histogram[ss2,{0.05},"PDF"],Plot[{Exp[-s],Pi s /2 Exp[-Pi s^2/4],32s^2/Pi^2 Exp[-4 s^2/Pi]},{s,0,4}],PlotRange->{{0,4},All}]

Show[Import["/home/broniows/kielce/Computer_Methods_2015/chaos/widma.png","PNG"],ImageSize->900]

Liczby losowe, liczby pierwsze, jądro złożone (podobnie dla widm atomowych, molekularnych), bilard Sinaia (GOE), rozkład dodatnich części urojonych biegunów funkcji ζ Riemanna (GUE), Laplasjan na rozmaitościach Riemanna, spektralna teoria grafów

Chaos algebraiczny

Quit[]

Ta dziedzina matematyki zajmuje się własnościami spektralnymi operatorów na przestrzeni Hilberta.

Przykładowo, zbadajmy rozkład odległości miedzy dalekimi liczbami pierwszymi sąsiednimi o di=1,2,3, ...

di=2

off=1000000
max=3000000
Prime[off]

Timing[tt=Table[Prime[off+i+di]-Prime[off+i],{i,1,max}];]

Mean[tt]//N

tt=tt/Mean[tt]//N;
Show[Histogram[tt,{0.21},"PDF"],Plot[{Exp[-s],Pi s /2 Exp[-Pi s^2/4],32s^2/Pi^2 Exp[-4 s^2/Pi]},{s,0,4}],PlotRange->{{0,4},All}]

Koniec! Dziekuję za uwagę - WB