Computer methods of modelling physical systems

Lectures for PhD studies, Institute of Physics, Jan Kochanowski U., Kielce

2017/18

Wojciech Broniowski

LOOK ONLY AT LECTURES MARKED (*)

Newest version of this file to be found at http://www.ujk.edu.pl/~broniows

Ver. 17.3, 17 April 2018, Mathematica 11.3

Setting the data path, modify for your case

(* Windows, example *)
(* datapath="d:\\my_dir\\"; *)

(* unix *)
datapath="/home/broniows/Dropbox/Computer_Methods_2018/"

$Path=AppendTo[$Path,datapath];

Global options for plots:

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

Lecture 1 (*)
Pseudo-random numbers and Monte Carlo integration

Built-in pseudo-random number generator

Properties of a good generator:

* lack of correlations (each number independent of the others)
* long period (before the sequence starts repeating)
* independemce of the choice of the seed
* speed

MLCG - Multiplicative Linear Congruential Generators

(Subscript[x, 0])=seed   - e.g. based on date and time
(Subscript[x, n+1])=((Subscript[x, n])*a+b) mod c  - a, b, c are properly chosen parameters

r=(Subscript[x, n])/c yields a number from [0,1],   c is very large, of the order of the period

For instance

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

r:=Random[]

Due to colon, at each call a diferent number is produced!

r

nsamp=10000
t=Table[r,{nsamp}];

The disthisribution is uniform (up to flucuations for the finite size of the sample).

nbin=100
his=Histogram[t,nbin,"PDF",FrameLabel->{"x","f(x)"}]

Show[his,Plot[{1+1/Sqrt[nsamp/nbin],1-1/Sqrt[nsamp/nbin]},{x,0,1},PlotStyle->{None},Filling->1->{2}]]

Monte Carlo integration

Function f(x), positive, set on the real range [a,b], bounded from above by M

f[x_]=Sin[x^2]^4
M=1
a=0
b=Sqrt[Pi]
Plot[f[x],{x,a,b},FrameLabel->{"x","f(x)"}]

Let us try to integrate analytically (!)

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

It works, and the result contains a special function (hitting F1 key will display more info on the Fresnel function). 
Let us use the built-in numerical integration:

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

Next, we will evaluate the same by and explicit Monte Carlo method. 
First, we make a generator of a point from the rectangle [a,b]x[0,M]
(do not forget a colon!)

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

po

If the point is below the curve, it is accepted (nin++), otherwise it is rejected. We make np samplings. 
The vaue of the integral is approximated as 

nin/np  M (b-a) 

(area of the rectangle time the fraction of success, i.e., being below the curve)

exact//N

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

Compared with the earlier exact result, not bad! In Monte Carlo, the relative error scales as 1/np, hence getting high accuracy is costly.

Abs[app-exa]/exa

Now a more “visual” variant, which stores the points

np=100

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},FrameLabel->{"x","f(x)"}];

Show[plof,plo]

Buffon’s needle (1733, solved 1777) - example of Monte Carlo and graphics

Quit[]

A famous example of “analog” throwing of a needle of length L on ruled paper with lines separated by d. Counting the ratio of cases where the needle crosses a line to all cases allows to approximate π (!)

First, analytic considerations:

L=1;
d=1;

The position of one tip of the needle along the axis perpendicular to the lines is x+L/2 cos(ϕ) (the other tip is at x-L/2 cos(ϕ))

Plot3D[x+L/2 Cos[ϕ],{x,-d/2,d/2},{ϕ,-Pi,Pi},AxesLabel->{"(Subscript[x, center])","ϕ","(Subscript[x, tip])"}]

The regions in x and ϕ where we have a crossing, x+L/2 Cos[ϕ]>d/2 or -d/2>x+L/2 Cos[ϕ], are visible in the plot below.

Plot3D[{-d/2,x+L/2 Cos[ϕ],d/2},{x,-d/2,d/2},{ϕ,-Pi,Pi},AxesLabel->{"(Subscript[x, center])","ϕ","(Subscript[x, tip])"},Mesh->True,PlotStyle->{{Opacity[.7],Gray},{Opacity[.7],Red},{Opacity[.7],Gray}}]

From symmetry, it is enough to consider the upper right corner

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

Next, we take any L, and d = 1 may remain as the scale

L
Clear[L,x]
L

Solve[x+L/2 c == d/2,x]  (* condition for the boundary in x, notation c=Cos[phi] *)

We need to evaluate ((∫_0)^π)dϕ ((∫_0)^((1-L cos ϕ)/2))dx =(∫_(-1)^1)((d cosϕ)/((d cosϕ)/(dϕ))) ((∫_0)^((1-L cos ϕ)/2))dx=(∫_(-1)^1)((d c)/(Sqrt[1-(c^2)])) ((∫_0)^((1-L c)/2))dx =i2  

and ((∫_0)^π)dϕ ((∫_0)^(1/2))dx = i20

Abs[D[Cos[phi],phi]] (* Jacobian, we will integrate over c instead of ϕ - see the following *)
jac=1/Sqrt[1-c^2]

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

i2=Integrate[jac i1,{c,0,1}]

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

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

The following ratio gives π:

2L/d i20/i2

i.e., for any d and L we have  π = 2 L/d x(number of all throws)/(number of throws with line cuts)

Now the same with Monte Carlo:

Quit[]
glob

L=1
d=1

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

Needle at angle phi and the center at (x,y) on the plane:

needle:={phi,x,y}

needle

Condition of cutting (geometry, as above)

cuts[ne_]:=Module[{m},m=Mod[ne[[3]],d];m=Min[m,1-m];L/2 Cos[ne[[1]]]>m]

cuts[needle]

gr[ne_]:={If[cuts[ne],Red,Blue],Line[{{ne[[2]]+L/2 Sin[ne[[1]]],ne[[3]]+L/2 Cos[ne[[1]]]},{ne[[2]]-L/2 Sin[ne[[1]]],ne[[3]]-L/2 Cos[ne[[1]]]}}]}

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

needles={} (* no needles yet *)

Do[needles=AppendTo[needles,gr[needle]],{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[needles],Axes->False]

needles={}

ListAnimate[Table[Show[canvas,Graphics[needles=AppendTo[needles,gr[needle]]],Axes->False],{i,1,300}]]

... finally, computing π in n throws

s=0 (* number of throws with cuts *)
n=1000000 (* all throws *)
Timing[Do[If[cuts[needle],s++],{i,1,n}];]
res = 2 L n/(d s)//N
100 Abs[(res-Pi)/Pi] "%"

Built-in distributions

Quit[]
glob

RandomVariate[NormalDistribution[0,4]]

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

Any distribution via inverse cumulative distribution function (CDF) - not done yet

Consider the transformation of probability distributions with the change of variabl
p(y) dy  =f(x)dx    
(but f(x)=1, as we assume it is uniform on [0,1]).  Hence  p(y) dy  =dx and CDF for y is
P(y)=∫p(y)dy =x
and
y=(P^(-1))(x) 
This simply meads that to obtain any distribution, one needs first to generate a uniform variable x and next plug it into the inverse CDF of the requested distribution

Example :

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,FrameLabel->{"x","f[x]"}]

Any distribution via the von Neumann method - not done yet

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->{Red,Thickness[0.01]}];
Show[plh,plf,FrameLabel->{"x","f[x]"}]

HOMEWORK

0) Install Mathematica and carry out tutorial. Run the lecture file.
1) Program in any multiplicative random number generator. 
2) Construct by a) inverse CDF, and b) von Neumann method the generator for the radial Woods-Saxon distribution f[r]=r^2/(1+Exp[r]). Compare, which method is faster. Hint: in part a), for speed, one needs to tabulate the inverse CDF and construct its interpolation via “Interpolation”.

Solution 1

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,{100000}];
Histogram[t,10,"PDF"]

Solution 2 a)

Clear[a]

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 (* normalized function *)

pl=Plot[fn[r]/.{R0->1,a->.1},{r,0,2},FrameLabel->{"r","4((πr)^2)WS(r)"}]

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

Plot[{F[y,1,.1],1},{y,0,2},FrameLabel->{"r","P(r) - CDF"}]

tt0=Table[{x,F[x,1,.1]},{x,0,2,2/1000.}]; (* table {x,F[x]} *)

tt=Table[{tt0[[i,2]],tt0[[i,1]]},{i,1,1000}];  (* table {F[x],x} *)

Finvint=Interpolation[tt]

Plot[Finvint[x],{x,0,.999},FrameLabel->{"y","(P^(-1))(y)"}]

w:=Finvint[Random[]] (* our basic formula *)
Timing[tab=Table[w,{100000}];]

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

Lecture 2  (Present 5) (*)
Ising model and the Metropolis algorithm

Introduction

Mini-primer of statistical physics:

A system is assumed to be in thermal equilibrium with the heat bath, h - external magnetic field, M - magnetization

Occupation probability of state i of energy (Subscript[E, i]) and magnetization (Subscript[M, i]): (Subscript[P, i])=1/Z (e^(-((Subscript[E, i])-(Subscript[M, i]) h)/((Subscript[k, B])T)))
Normalization (partition function): Z(T,h)= (∑_i)(e^(-((Subscript[E, i])-(Subscript[M, i]) h)/((Subscript[k, B])T))) - contains all statistical information on the system,
β=1/T, the Boltzmann constant (Subscript[k, B])=1

(Kerson Huang)

mean energy:<E>=(∑_i) (Subscript[P, i]) (Subscript[E, i])= - (∂/(∂β ))log Z - (∂/(∂β ))log (∑_i)(e^(-((Subscript[E, i])-(Subscript[M, i]) h)β))

specific heat: (Subscript[C, V]) ≡ ((∂<E>)/(∂T)) = (β^2) [(1/Z) (((∂^2)Z)/(∂(β^2))) -  (1/(Z^2)) (((((∂Z)/(∂β))))^2) ]= (β^2) [<(E^2)> - <E((>)^2)] = β^2((ΔE)^2) - fluctuation-dissipation theorem (!)

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

magnetic susceptibility: χ ≡ ((∂<M>)/(∂h))=β [<(M^2)> - <M((>)^2)] 


Magnetics as spin networks

* A simple model of magnetism 
* Example of MC simulation in statistical physics
* Model of phase transitions/critical phenomena, extension to quantum field theory (lattice-gauge calculations) 

Magnetics undergo phase transitions. Magnetic susceptibility becomes divergent at the Curie temperature

	low temperature ⟶ order ⟶ a common direction of spins  ⟶ high magnetization
	high temperature ⟶ disorder ⟶ random directions of spins ⟶ no magnetization

1D Ising model (Wilhelm Lenz 1920, Ernst Ising 1925 )

Energy of the system: E=-J  (∑_(<i,j>))(Subscript[S, i])(Subscript[S, j])
(Subscript[S, i])= +1 lub -1   - spin at node i of the lattice 
<i,j>  -  closest neighbors
J   - interaction strength

When J > 0  - ferromagnetic, if J < 0 - antiferromagnetik (we will consider the ferromagnetic case)

Ising: no phase transition in 1D
Lars Onsager (1944): exact solution of the 2D case, existence of the phase transition as in real ferromagnetics
3D case - not solvable analitically, simulations

Back to statistical physics

Magnetization is a thermodynamic average over states C, 

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

T - temperature
C - configuration
P(C,T) - probability of configuration C
M(C) - magnetization of konfiguration C

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

p(C)=(Z^(-1)) (e^(-E(C)/T))
Z=(∑_C) (e^(-E(C)/T)) - partition function
E(C) - energy of configuration C

Monte Carlo method 

Naive method based on independent sampling

<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))

will not work! Most of the randomly chosen configurations will have negligible weights, there is no chance of hitting “good” configurations

Brute-force method does not work

(Subscript[N, A ])=6(x10^23) - the Avogadro number of spins in a macroscopic systems. A typical 2D lattice has much less: 24x24, hence the number of configurations is  (2^(24x24)) ≈ (10^172)
(10^(9 ))configurations/processor/s x(100^( ))processors x (10^14)s/year x 1year  = (10^25) - we can check only this many (a 9x9 lattice)!

Importance sampling

If we sampled configurations having the Boltzmann distribution (as in the thermodynamic sum), then the Monte Carlo mean would be

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

Unfortunately, we do not have the distribution - we are trying to compute it! 

Markov chains

P(A,t) - probability of configuration A after “time” (number of steps) t
W(A→B) - transition probability per time unit from A to B 

The gain-loss equation:

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

At large t we have 

P(A,t)→p(A)     (**)

- loss of memory of the initial conditions

Detailed balance
 
With (**), the sufficient condition for (*) is 

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

or

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

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

The Metropolis algorithm (Metropolis, Teller, Rosenbluth, 1953)

How to satisfy (***)? One (of infinitely many) possibilities:

W(A→B)=(e^(-(E(B)-E(A))/T))<1       if  E(B)>E(A)  (excitation with probability exponentially decreasing with the excitation energy)
                   =1                           if  E(B)⩽E(A)  (always goes to a state with a lower energy)

Then if E(B)>E(A), 

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

whereas for E(B)<E(A) 

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

- the same! (as should be)

General requirements

Monte Carlo algorithm must fufill

* the detailed balance - OK for Metropolis
* ergodicity (any configuration must be pssible to reach from another one in a finite number of steps) - depends on the model, temperature - generically, a highly non-trivial problem

The Metropolis algorithm for the 1D Ising model

Quit[]

(* number of lattice nodes *)
n=100

(* coupling constant, +1 ferro, -1 ferri, we assume |J|=1 as a scale *)
J=-1

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

init

(* test of periodic "winding" - convenient below to write sums in a uniform way *)
Table[{i,Mod[i+1,n,1]},{i,1,n}]

(* spin, energy and spin correlation (per node) *)
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

(* m - displacement *)
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]}

(* operation of spin flip in node i *)
flip[i_]:=t[[i]]=-t[[i]]

(* change of energy at spin flip at node i *)
del[i_]:=2J t[[i]](t[[Mod[i-1, n,1]]]+t[[Mod[i+1, n,1]]])

(* The Metropolis jump in node i at temperature T *)
step[i_,T_]:=Module[{},If[del[i]<0,flip[i],If[Exp[-del[i]/T]>Random[],flip[i]]]]

(* a single "sweep", starting from a randomly selected nore 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}];]

Let us see how it works

init

energy

(* try a few times *)
sweepT[.1]
energy

(* reaching the thermodynamic equilibrium at temperature T, number of sweeps = 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}]*)]]

(* visualization of the state - try different temperatures *)
init;
Show[Graphics[Raster[{(t+1.2)/3}]],ImageSize->800,Background->LightBlue]
etrace[.2, 1000]
Show[Graphics[Raster[{(t+1.2)/3}]],ImageSize->800,Background->LightBlue]

(Legended[, (PointLegend["energy/n", "spin/n"])])

Simulations

ninit=200; (* number of sweeps needed to achieve equilibrium ... *)
ncoll=1000; (* ... and the number of following "measurement" sweeps  *)

(* average (over ncoll measurement states) of spi, energy, specific heat, and magnetic susceptibility *)

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]

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

(* plots *)
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]}]

Show[ListPlot[tCV,PlotStyle->Red],ListPlot[tchi],ImageSize->600,AxesLabel->{StyleForm["T",FontSize->18],StyleForm["(Subscript[C, V])", "χ",FontSize->18]},PlotRange->All]

Exact Ising’s solution

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

(* comparison of Monte Carlo and the exact cases *)
Show[plmc,Plot[ene[T,n],{T,0,5}]]
Show[plmc2,Plot[CV[T,n],{T,0,5}]]

Finer grid in T for (Subscript[C, V]) and χ

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}]]

(* exact result for χ at large T *)
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}]]

Correlations

(can also be obtained analytically with the help of the so-called transition matrix)

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]

(* from construction, gamma[0]=1 *)
gamma[0]

(* matching to an exponential function *)
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]]

HOMEWORK

1) Redo this chapter for the ferrimagnetic case. What are the differences compared to ferro?
2) Read the notes by Wei Cai (Caltech), file Chap12_Ising_Model_v04.pdf
3) Program the Metropolis simulation for (choose one case)

 - 1D Ising with h≠0
 - 1D Ising with 2 or more rows
 - 1D Ising with interactions farther than with the nearest neighbor
 - 1D Potts model (a generalization of the Ising model for (Subscript[s, i])=1,2,..,q), (Subscript[H, ij])=-J cos(2π((Subscript[s, i])-(Subscript[s, j]))/q)
 - 2D Ising
 - The (classical) Heisenberg model in 3D, (Subscript[H, ij])=-J (OverVector[(Subscript[s, i])])·(OverVector[(Subscript[s, j])]) ((OverVector[(Subscript[s, i])]) is a classical vector)
 
The template for above is the present chapter, which may be modified appropriately

Lecture 3 (present 6)(*)
Further applications of the Metropolis algorithm: 
from combinatorics to path integrals

Design of integrated circuits

Consider the following problem: we have n transistors. The designed logic of the system determines the connections between them, i.e., which pairs are connected. The capacity (size) of the silicon plate is limited; suppose we need to use two plates to build the system, as we need more transistors that can fit on one plate. The transistor connections within the same plate are cheap, whereas connections between different plates are expensive (as they require wires). Goal: divide the transistors betwen the two plates in such a way that the number of cross-plate connevtions is minimized. 

Solution: this is a hidden spin problem. Assign to transistor i label s(i)=0 if it belongs to the first sheet, and s(i)=1 if it belongs to the second sheet. The “Hamiltonian”  (cost function) can be written as 

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

Next, we proceed according to the annealing algorithm, i.e., after each Metropolis sweep we lower the temperature a bit.

Simulation

Quit[]

(* number of transistors *) 
n=50
(* maximum number of connections *)
k=n(n-1)/2

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

rs

(* number of possibilities = n over n/2 *)
n!/(n/2)!/(n/2)!//N

Binomial[n,n/2] //N

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

init
t

(* randomly generated connections *)
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}]

(* swap of transistors i and j *)
flip[i_,j_]:=Module[{},to=t[[i]];t[[i]]=t[[j]];t[[j]]=to;]

(* the Metropolis jump *)
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

In our annealing procedure the temperature decreases linarly at each sweep form the value T down to T/(total number of sweeps). Of cource, the desired limit T ⟶ 0 may be obtained in another way.

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
energy
etrace[.35,100]
t
energy

Annealing, frustration, and spin glass

(* on blackboard *)

Feynman path integrals

From quantum mechanics we recall that the evolution operator can be represented as an integral over trajectories (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]))/ℏ)),

where A is a constant and S=(∫_(Subscript[t, i])^(Subscript[t, f])) dt L[x(t),(OverDot[x])(t)] denotes the classical action for a particle with trajectory x(t), with fixed x((Subscript[t, i]))=(Subscript[x, i ])oraz x((Subscript[t, f]))=(Subscript[x, f ]). We assume the units ℏ=1. Discretizing over time, with the velocity v=((x((Subscript[t, j+1]))-x((Subscript[t, j])))/(Δt)), we get

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]))])],

where (Subscript[t, 1])=(Subscript[t, i]), (Subscript[t, n])=(Subscript[t, f]), Δt = ((Subscript[t, f])-(Subscript[t, i]))/n. Next, we pass to imaginary time 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])])],

with Τ=(Subscript[τ, f])-(Subscript[τ, i]), Δτ=T/n, (Subscript[x, j])=x((Subscript[t, j])). Formally, the above expression is analogous to a spin model in 1D, where the value of “spin” is replaced with a coordinate x (real number). The interaction between the neighboring sites comes form the kinetic energy, and the inverse temperature is proportional to the Euclidean time,  β∼Τ. 

The normalization (for a free case A = (((2 ((πΤ)^2)/n))^(-n/2))) does not depend on (Subscript[x, i]), hence is irrelevant.

Below we consider a case of a particle of mass m=1 in the harmonic potential V(x)=(x^2)/2.

Generally, QM (0 space dimensions, 1 time dimension, where we vary the function x(t), is equivalent to a 0+1 D spin model, whereas a field theory in d+1 dimensions (the variable is a field ϕ(t,x)) is equivalent to a D=d+1 -dimensional spin model (lattice QCD etc.)

PIMC :  Path-Integral Monte Carlo

Quit[]

(* number of latice sites *)
n=50

(* limit for the coordinate *)
sc=3;

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

(* random initialization of the trajectory *)
init:=Module[{},t=Table[rn,{i,1,n}];t[[n]]=t[[1]];]

init
t

(* energy E_kin+E_pot *)

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

(* energy from the virial theorem - for the harmonic oscillstor E_kin=V, E_tot=2V - this yields better stability *)
ve[tau_]:= 2Sum[t[[i]]^2/2*(tau/n),{i,1,n}]/tau

(* operation of a small change of the coordinate in site i *)
flip[i_]:=t[[i]]+=rn/10

(* the Metropolis jump in site i at the Euclidean time τ *)
(* we impose periodicity 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]];]

(* a sweep consisting of one jump only *)
sweepT[tau_]:=Module[{r},r=Random[Integer,{1,n-1}];step[r,tau];]

Illustration of trajectory smoothing

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]

PIMC simulations

We will find the ground-state energy and  the modulus square of the wave function for the  harmonic oscillator (HO).  Of course, for HO
(Subscript[E, 0])=1/2 and |(Subscript[ψ, 0])((|)^2)=(e^(-(x^2)))/√π.  On the other hand, in the energy basis
 
 Z=Tr U(T)=
∫dx<x| U(Τ)|x>
=∫dx (∑_n)<x| n > < n | U(Τ) | n > < n |x>,
where
< n | U(Τ) | n > = (e^(-(Subscript[TE, n]))) 
(note that in QM we have  U=(e^(-TH)), where H is the time-independent Hamiltonian of the system). 
At large Τ we have (excited states are  "killed")

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

Hence (Subscript[E, 0])=(lim_(Τ →∞))Log[-Tr U(T)]/Τ, i.e., the free energy  tends to the energy of the ground state.  

Another method (used in simulations below): 

The mean for the time-independent operator O at fixed time Τ is 

<O(Subscript[>, Τ]) = ((Tr [O U(Τ)])/Z) = (((∑_n)<n|O|n> (e^(-Τ (Subscript[E, n]))))/((∑_n)(e^(-Τ (Subscript[E, n])))))⟶ <0|O|0> for large Τ,
i.e., it approaches the ground state expectation value. 

Beacuase in the Metropolis method the states are probed with the  Boltzmann-like distribution (e^(-S)),where S denotes the action for a given trajectory, in the example below the sums are carried out weigh 1. Hence, at large Τ for  O=H we get the ground-state energy,and for O=|x><x| the modulus squared of the graound-state wave function. 
Other observables may be evaluated in a similar way.

ninit=5000; (* number of (rejected) sweeps until the equilibrium is reached ... *)
ncoll=100000; (* ... and the number of the following measurement sweeps *)

(* average (over ncoll trajectories) energy from the virial theorem *)
eav2[tau_]:=Module[{},init;Do[sweepT[tau],{k,1,ninit}];scoll=0;Do[sweepT[tau];scoll+=ve[tau],{k,1,ncoll}];scoll/ncoll]

(* average (over ncoll trajectories) modulus squared of the wave function *)
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.]

(* symmertization for better statistics *)
psi=Join[psi,-psi];

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

Mean[psi]
(* variance should be 1/2 for the ground state *)
vp=Variance[psi]

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

(* energy of the ground state *)
e0=1/2

eav2[15.]

(* energy as a function of τ (=T) *)
Timing[teav=Table[{tau,eav2[tau]},{tau,.5,50,3}]]

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}]

No improvement from increasing τ beyond a certain value!

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

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

The time resolution a=Τ/n must smaller than a characterisctic scale, a=T/n<<1/(Subscript[E, 0]), whereas
Τ must be larger than the inverse energy gap (to suppress excited states), Τ= n a>> 1/((Subscript[E, 1])-(Subscript[E, 0])).

Obtaining information on excited states is costly.

Improved actions are such discrete actions for which the continuum limit a ⟶ 0 is faster (i.e., achieved at lower values of n).

Lattice gauge theory: vertices - coordinates (t, (x^1),(x^2), (x^3)), variables in vertices - fermion fields, links - gauge fields.

HOMEWORK

1) Read “Computational Quantum Physics” by de Forcrand and Troyer, chapters 5, 7-9.
2) Modify and run the code from the lecture for another potential.

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 β

Lecture 6 (present 2) (*)
Neural networks I

Instead of introduction

Mach’s illusion

Quit[]
glob

(* intensity function - gradual step *)
low=0.2; up=0.8;
xl=300; xu=600;

(* definition of function with conditions *)
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["intensity",FontSize->18]}]

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

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

One can note a lighter and a darker stripe. This is clearly an illusion, as the intensity function is monotonic! This proves some “image correction” algorithm working in our vision system.

Hermann’s illusion

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]

(Subscript[(Subscript[\[Placeholder], \[Placeholder]]), \[Placeholder]])

Leaf wave

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

Rotating snakes

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

Remarks

Our visual system (retina, brain) carries out image correction whicg leads to side effects.
The effect is dynamic (waving), the time scale is unexpectedly long, about 1s.
It applies to everyone, hence occurs od deep level  (wetware).

Neural networks - basic information

Show[GraphicsGrid[{{Show[Import["/home/broniows/Dropbox/Computer_Methods_2018/neuron/neuron.png","PNG"],ImageSize->300],
Show[Import["/home/broniows/Dropbox/Computer_Methods_2018/neuron/neural_network.png","PNG"],ImageSize->200]}}]]

The simplest “forward” neural network: input layer (dendrites), hidden layer(s) - can be many, output layer (axons). There is a very large number of connections between the neighboring layers (axon-dendrite).

The human brain contains about (10^11) neurons, where each has on the average about 7000 connections to other neurons. A single neuron works in a very simple manner, as a threshold transceiver, however the behavior of the whole brain is very complex.

Emergence of “complexity from simplicity”.

A simple model of the neuron (Warren McCulloch, Walter Pitts, 1942)

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

(Subscript[x, i]) - input signals, (Subscript[w, i]) ∈ℝ - weights, (Subscript[w, 0 ])- bias, f - activation function, (Subscript[w, i])>0 - excitation, (Subscript[w, i])<0 - damping

f[S_]=UnitStep[S];
Plot[f[S],{S,-2,2},PlotStyle->Thickness[.01],AxesLabel->{StyleForm["signal",FontSize->18],StyleForm["does not fire / fires            ",FontSize->18]},ImageSize->600]

Let us note that the activation function f(x) is nonlinear, i.e., f(x+y) ≠ f(x)+f(y), f(a x)≠a f(x). Nonlinearity is a key property of neural systems. Because nonlinearity is not physical, one frequently uses as smooth function, the so-called“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]

MCP neuron (with a sharp threshold):

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

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

Logical gates and networks

With MCP neuron we can form logical networks

(* AND gate *)
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}];

(* NAND gate *)
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}];

(* OR gate *)
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/Dropbox/Computer_Methods_2018/neuron/xor.jpg","JPG"],ImageSize->400],Pi/2]

(* we move backward to subsequent layers of the network *)
neurXOR[x1_,x2_]=neurAND[neurNAND[x1,x2],neurOR[x1,x2]];

(* XOR gate *)
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]

We can see that the basic logical gates can be formed. From these gates one can form any logical network (gates NAND and XOR are complete, in the sense that any network can be formed only with them). A multilayer MCP network is thus equivalent to a universal Turing machine.

Pattern recognition

Show[Import["/home/broniows/Dropbox/Computer_Methods_2018/neuron/TH.png","PNG"],ImageSize->900]

The Hamming rule - closest “distance” from the pattern

upper row                 middle row                         lower row

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

The columns at Xij correspond to the rows in the figure, first columns for T, second for H. The values Fi are the Hamming functions.

Digression

aa={1,2,3}

aa//FullForm

aa[[0]]

ff=Apply[g,aa]

ff[[0]]

d[s1_,s2_]:=Apply[Plus,(s1-s2)^2]  (* distance *)

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}]]

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

We divide the answers into classes based on the distance form the pattern:

recognize[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,"this is T",test<3 && test >=2,"rather T", test>1 && test<2,"do not know",test≤1 && test >0,"rather H",test==0,"this is H"]]]

recognize[{1,1,1},{0,1,0},{0,1,0}]
recognize[{1,0,1},{1,1,1},{1,0,1}]
recognize[{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]

recognize[xr1,xr2,xr3]

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

Linear neuron as  a memory model

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}};

something=
      {{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[something]}},ImageSize->800]

We unfold the above matrices into vectors for an easier algebraic processing:

fA=Flatten[A];
fa=Flatten[a];
fI=Flatten[II];
fi=Flatten[ii];
fsomething=Flatten[something];

Heteroassociative memory

Digression from linear algebra:
The outer (diadic) multiplication of vectors makes a matrix:

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

We have the desired projection features 

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

if w.w’=0, i.e., if vectors w and w’ are orthogonal,

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

For a large dimensionality n, randomly selected k vectors, k<<n, will be approximately orhogonal. This feature is used in the construction of pattern recognition explained below.

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

Now we make a “linear model of heteroassociative” memory.
(“Partition” is a reverse of Flatten, making a matrix from a vector)

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}]

We get an association between pairs a-A, i-I. Small errors ollow from a not exact orthogonality.

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

Let us try to get rid of the errors be applying a suitable filter which will remove the lighest shades. Floor[x] yields the highest integer less or equal to x. In addition, we apply “Unitize”, which replaces all integers >0 with 1. For instance, for “i” we get:

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

This (nonlinear!) filter evidently helps! Thus we have created a good model of heteroassociative memory.

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

Autoassociative memory

Next, we associate images with themselves (autoassociation)

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]

Undistorted image before and after filtering:

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]

Distorted image, i.e.,  a number “bad” of bad bits is randomly generated. The recognition still works good if “bad” is not to large. The filter threshold is relevant.

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

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

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

We get the error correction from a distorted/incoplete input!

guess[ba_]:=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}],{ba}];
GraphicsGrid[{{ArrayPlot[Partition[fAcor,12]] (*,pl[fAcor],plf[po fAcor]*)}},ImageSize->200]]

guess[30]

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"]

Digression on linear networks

We have just seen that linear networks (in particular, after application at the end of a nonlinear filter) can do useful calculations. The structure of such networks in the neural interpretations looks as follows: they have two layers, input i output, whose neurons are connected with edges of weights (Subscript[w, ij]). Their functioning is equivalent to matrix multiplication

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

Because of linearity it does not make sense to consider networks with more layers, as then we would have

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

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

Hence linear multilayer networks (with a possible nonlinear filter, e.g. a “sigmoid”, applied at the end) are equivalent to a two-layer linear network (input-output + possible filter). Such networks cannot do all computations, in particular, they cannot implement the XOR gate.

Part of neocortex responsible for vision

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

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

Homerk

1) Read the paper M. P. van den Heuvel, O. Sporns, “Rich-Club Organization of the Human Connectome”, The Journal of Neuroscience, 31(44):15775 (2011) - file 15775.full.pdf

Lecture 7 (present 3)(*)
Neural networks II

Supervised learning of the neural network

Perceptron (Frank Rosenblatt, 1957)

input layer            hidden layers                                         output layer

Simple neural network with input, hidden, and output layers.
NB, the appearance of the book “Perceptrons: an introduction to computational geometry” (Marvin Minsky, Seymour Papert, 1969) led to the “AI winter”. Multi-layer networks do not have thereby discussed limitations.

Two stages: 1) learning 2) using.

Let the network with weight {w} takes an input vector x and yields an output vector y, i.e., y=S({w},x). 
Let the pairs (x,t) denote input and the desired (correct) output (we must know the proper answers to train the network). There anr N such pairs and they form our training sample.
  
1) For each training sample we carry out the (supervised) learning, based an adjustment of the weights that the error in the output is minimized. We do it, e.g., as follows: in a loop over the inputs in the sample, a=1..N, we take (Subscript[x, a]) and propagate it forward, obtaining output (Subscript[y, a])=S({w},(Subscript[x, a])). We evaluate the error, usually in the form 

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

The error depends, of course, on the selection of the weights. We thus treat the error (for a given training sample) as a fuction of {w} and simply minimize E({w}). However, this may not be a simple task because of the presence of local minima.

A popular algorithm is the presented below backward propagation, based on the the steepest descent. It updates the weights ater each training step, going backwards from the last to the first layer. 
 
2) Having set the weights (trained the network), we have a network dedicated to a given task. For a given x we get the answer y=S({(Subscript[w, min])},x). 

Remarks:
 
The perceptron a certain nonlinear interpolation. 
“Programing” is equivalent to finding appropriate weights.
Local minima are apparently not a problem for large multi-layer networks (deep learning).
More sophisticated methods of training: annealing, genetic algorithm...)
Larger networks have more power, but need a longer training.

Sigmoid

Quit[]

(* sigmoid *)
f[z_]=1/(1+Exp[-z]) 

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

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

For the sigmoid f’(x)=f(1-f)

a1==a2 //Simplify

Sometimes one introduces the sigmoid with a “temperature”, 
f(x;T)=1/(1+Exp[-x/T]). Parameter T washes out the  “decision range”. In our examples we set  T=1.

Backward propagation

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

We now construct a perceptron with 2 hidden layers

ni=3; (* number of input neurons *)
n1=7;  (* number of neurons in the first hidden layer *)
n2=7;  (* number of neurons in the second hidden layer *)
no=1;   (* a single output neuron *)

Conditions for three considered examples - read-in only one (!)

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

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

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

(* This case falls into a local minimum and the network does not work! *)
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

Random initial weights from the range (-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}];

W1//MatrixForm
W2//MatrixForm
W3//MatrixForm

Backward error propagation algorithm (Arthur E. Bryson, Yu-Chi Ho, 1969). From mid 70’s introduced to neural networks, which led to AI renaissance.

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=20000; (* number of learning steps must be sufficiently large *)

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

eta=0.1; (* learning "speed" *)

y1=f[x.W1]; (* flow vectors y in subsequent layers *)
y2=f[y1.W2];
y3=f[y2.W3];

(* Bryson-Ho backward propagation algorithm *)

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}]

W3 //MatrixForm

Now use the trained network ...

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

(Return to this cell later and repeat)

(* cutting some connections *)
W1[[3,4]]=0;
W1[[3,7]]=0;
W1[[3,6]]=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],{6000}],gr,ImageSize->700]

Quiz

How do we manage “by naked eye”?

x0={Random[],Random[],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]

Brain damage

Graphic visualization

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]

Now return higher

For other initial conditions and training sequence the network can be different!

Summary

Remarks on the multi-layer perceptron:

- Long training on a large sample improves results.
- Same sample may be reused in iterations of the backpropagation algorithm.
- Large networks, with many layers and neurons, have larger possibilities.
- The trained networks “doing the same job” may have completeley different sets of weights (local minima), though work with similar efficiency. 
- Connections with small weight magnitudes may be cut with no signicant deterioration of functionality (brain damage).
- The problem of falling inlto local minima may be circumvented with more sophisticated methods.

Let us also note that the above implemented figure recognition has nothing to do with geometry per se. It is purely “neural”.

Why it works? 

A single MCP neuron at a sharp threshold yield the condition (Subscript[w, 0])+(Subscript[w, 1])(Subscript[x, 1])+..+(Subscript[w, n]) (Subscript[x, n]) = 0, which is a hyperplane equation. Hence one may divide the points into two cathegories: above and below the hiperplane. In general, a multi-neuron perceptron składa combines (many) such conditions, thus being able to recognize if a point belongs to a more complicated region.

HOMEWORK

0) Read “Neural Networks” (Christos Stergiou, Dimitrios Siganos),  file nn.pdf.
1) Verify if the network presented in the lecture works for other figures, also figures with holes, e.g. ring in a square. More layers may be needed.
2) Recompute the formulas of the backpropagation algorithm based on the file back_prop.pdf.

Lecture 8 (present 4)(*) 
Neural networks III

Crab’s eye / lateral damping
Hopfield networks

Remember Mach’s illusion?

(* intensity function *)
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/Dropbox/Computer_Methods_2018/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/Dropbox/Computer_Methods_2018/neuron/Hopfield-net.png","PNG"],ImageSize->300]

Model with feedback for N neurons in the same layer:

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

(Subscript[f, i]) - signal after i’th iteration (vector of dimension equal to the number of neuronów in the layer, i.e., N)
s - signal (vector of dimension N)
ϵ - paramter
W - symmetric coupling matix of dimension N x N

For a steady state (after many iterations) we have (Subscript[f, n]) ⟶ f, so
 
 s+W.f - f=0,
 f=(((1-W))^(-1)) s.         (***)

Let us take a sample form of the coupling matrix 

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

i.e., the dampling  tłumienie is strongest for close neighbors, and falls with distance with a characteristic skale b.

ns=30;       (* number of neurons, N *)
b=6;           (* damping fall-off parameter *)
a=0.2;     (* strength of damping *)
iter=100;  (* number of iterations *)
eps=.1;    (* "learning speed" parameter *)

(* a sample signal - Mach's paradox *)
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["input signal",FontSize->18]}]

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

(* neuron does not damp itself! *)
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["Coupling weights W",FontSize->18]]

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

(* implementation of a single iteration *)
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,FrameLabel->{"x","response in subsequent iterations"}]

rescale=2.5; (* for a better comparison *)
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.15}]],PlotRange->All,ImageSize->600,FrameLabel->{StyleForm["x",FontSize->18],StyleForm["input and obtained signals",FontSize->18]}]

One can see that in the crab’s eye model the flat fragments of the input signal are distorted. Qualitatively, this corresponds to Mach’s illusion. In particular, for x=10 we have a minimum (the lightest stripe), and at x=20 a (local) maximum (the dark stripe). Hence, a very simple model of negative feedback between the receptors in the retina is sufficient to explain the phenomenon qualitatively. 

Why does Natute take so much effort to make our senses deceiving? Answer: the same mechanism of image transformation leads to enhance vision contrast.

(* another sample signal *)
s=Table[3^2/((x-15.)^2+3^2),{x,1,30}];

g0=ListPlot[s,PlotRange->{{0,30},{-0.1,1.1}},ImageSize->400,FrameLabel->{StyleForm["x",FontSize->18],StyleForm["input signal",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.05}]],PlotRange->All,ImageSize->600,FrameLabel->{StyleForm["x",FontSize->18],StyleForm["input and obtained signals",FontSize->18]}]

The tansformed signal (solid line) is narrower (better contrast) from the original (dots), hence the network has carried out a useful image correction.

Comments

For Hopfield networks (with symmetric weights, (Subscript[w, ji])=(Subscript[w, ij ])) one can introduce the energy

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

Conditions for the minimum yield 

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

(moreover, matrix w must by positive-definite), from where f=(w^(-1))s (the same solution as (***) for w=1-W). Therefore, the Hopfield networks are equivalent to a spin model with couplings (Subscript[w, ij]) in an external field s! The quantity (Subscript[s, i]) may also be interpreted as the theshold of the MCP neuron.

Beacause the above condition has the matrix form f=M s (with M=(w^(-1))), it falls into the class of the earlier-discussed models of heteroasociative and autoassociative memory (the difference is only in the matrix inversion and the related iterations. Thus, the Hopfield networks serve also to model the neural memory.

Boltzmann machine (network)

This is a generalization of the symmetric Hopfield network for 

(Subscript[f, i])=0,1 (two states), 

which introduces a thermal (statistical) noise. Formally, the statistical partition function is 

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

where a denotes the configuration of the whole system and (Subscript[E, a]) is its energy, defined above (see also the lectures on the spin network,Metropolis, annealing, ...). One may also introduce visible and hidden nodes (neurons). The training of the Boltzmann machines is difficult for a larger number of neurons and strong coupling, which is related to thermal activation of neurons. A positive aspect is the fact that the temperature allows to go out of the local minima (as in annealing).

Homework 4

0) Read the first part of chapter 15 from the book by Raúl Rojas, Neural Networks, a Systematic Introduction (neuron_book.pdf).

1) Explain qualitatively the Hermann illusion, based on lateral dampling. What can you infer on the functioning of the yellow spot and the other receptors of the retina?

Lecture 9 (present 5) 
Neural networks IV

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

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

Wykład 10
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 11
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 12
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 13
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
b=1

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->600]

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

Eigenvalues[A]

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=1

off=1000000
max=4000000
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