Einführung in Mathematica XVIII

Metropolis-Algorithmus

© Jürgen Schnack, 2009/06/18

Problem

Approximation statistischer Mittel durch Mittelung über 
geeignete Mengen von Mikrozuständen. Erzeugung dieser Mengen durch den Metropolis-Algorithmus.

Beispiel: Ising-Modell

Wechselwirkungsmatrix

NMax=10;
NMaxSq=(NMax^2);

Spins sind wie folgt abgezählt:
(({{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}))

Hilfsfeld, in dem die Spinnummer in den Platz umgerechnet wird

Platz=Table[{IntegerPart[((k-1)/NMax)]+1,Mod[k-1,NMax]+1},{k,1,NMaxSq}]

Wechselwirkungsmatrix JMat

JMat=Table[0,{i,1,NMaxSq},{j,1,NMaxSq}];
For[i1=1,i1<NMax+1,i1++,
For[j1=1,j1<NMax+1,j1++,
For[i2=1,i2<NMax+1,i2++,
For[j2=1,j2<NMax+1,j2++,
(* Print[i1," ",i2," ",j1," ",j2," ",(i1-1)*NMax+j1," ",(i2-1)*NMax+j2]; *)
If[Mod[Abs[i1-i2],NMax]==1&&j1==j2,JMat[[(i1-1)*NMax+j1,(i2-1)*NMax+j2]]=1];
If[Mod[Abs[i1-i2],NMax]==NMax-1&&j1==j2,JMat[[(i1-1)*NMax+j1,(i2-1)*NMax+j2]]=1];
If[Mod[Abs[j1-j2],NMax]==1&&i1==i2,JMat[[(i1-1)*NMax+j1,(i2-1)*NMax+j2]]=1];
If[Mod[Abs[j1-j2],NMax]==NMax-1&&i1==i2,JMat[[(i1-1)*NMax+j1,(i2-1)*NMax+j2]]=1]
]
]
]
];
MatrixForm[JMat]

Ausgangszustand

Spins=Table[2 RandomInteger[]-1,{i,1,NMax},{j,1,NMax}];
Spins0=Spins;
Spins

I1=Show[Graphics[Raster[Spins]/.{1->{1,0,0},-1->{0,0,1}},AspectRatio->1]]

Metropolis für eine Temperatur

MaxSweeps=100;
T=2.6;

Im Verfahren wird Spin Nr. i umgeklappt (single spin flip).
DeltaE=(ENeu-EAlt) ist zweimal die Summe, weil man -s_i*Summe-s_i*Summe hat.

Spins=Spins0;
For[Sweeps=1,Sweeps<MaxSweeps+1,Sweeps++,For[i=1,i<NMaxSq+1,i++,DeltaE=2 Spins[[Platz[[i,1]],Platz[[i,2]]]] (∑_(k=1)^NMaxSq)JMat[[i,k]] Spins[[Platz[[k,1]],Platz[[k,2]]]];If[DeltaE≤0,Spins[[Platz[[i,1]],Platz[[i,2]]]]=-Spins[[Platz[[i,1]],Platz[[i,2]]]],If[RandomReal[]≤(E^(-(DeltaE/T))),Spins[[Platz[[i,1]],Platz[[i,2]]]]=-Spins[[Platz[[i,1]],Platz[[i,2]]]]]]];
Print[Graphics[Raster[Spins]/.{1->{1,0,0},-1->{0,0,1}},AspectRatio->1]]];

Neuer Ausgangszustand

Spins=Table[0*2 RandomInteger[]-1,{i,1,NMax},{j,1,NMax}];
Spins0=Spins;
Spins

Graphics[Raster[Spins]/.{1->{1,0,0},-1->{0,0,1}},AspectRatio->1]

Metropolis für verschiedene Temperaturen

MaxSweeps=100;
TSteps=20;
Magnetisierung=Table[{0.0,0.0},{k,1,TSteps}];

DeltaE=(ENeu-EAlt) ist zweimal die Summe, weil man -s_i*Summe-s_i*Summe hat.

For[iT=1,iT<21,iT++,T=0.5` iT;Spins=Spins0;Magnetisierung[[iT,1]]=T;For[Sweeps=1,Sweeps<MaxSweeps+1,Sweeps++,For[i=1,i<NMaxSq+1,i++,DeltaE=2 Spins[[Platz[[i,1]],Platz[[i,2]]]] (∑_(k=1)^NMaxSq)JMat[[i,k]] Spins[[Platz[[k,1]],Platz[[k,2]]]];If[DeltaE≤0,Spins[[Platz[[i,1]],Platz[[i,2]]]]=-Spins[[Platz[[i,1]],Platz[[i,2]]]],If[RandomReal[]≤(E^(-(DeltaE/T))),Spins[[Platz[[i,1]],Platz[[i,2]]]]=-Spins[[Platz[[i,1]],Platz[[i,2]]]]]];Magnetisierung[[iT,2]]=Magnetisierung[[iT,2]]+Abs[(∑_(k=1)^NMaxSq)Spins[[Platz[[k,1]],Platz[[k,2]]]]]]];Magnetisierung[[iT,2]]=((Magnetisierung[[iT,2]])/(Sweeps NMaxSq NMaxSq));Print[Magnetisierung[[iT,1]],"  ",Magnetisierung[[iT,2]]]];

ListPlot[Magnetisierung,PlotRange->{{0,10.5},{0,1.1}},PlotStyle->{Hue[0],PointSize[0.03`]},BaseStyle->{FontFamily->"Helvetica",FontSize->18}]

Schluß

Hier ist Schluß.