AB_MC.html

> restart;
with(plottools):
with(plots):

Warning, the names arrow and changecoords have been redefined

> n:=10;
omega:=-0.5;
Temp:=1;
AB:=Array(1..n,1..n);

> toss:=rand(0..1):

> for i from 1 to n do
for j from 1 to n do
AB[i,j]:=toss();
end do;
end do;

> AB;

> AB_display:=proc(AB)
local red,blue,i,j,pred,pblue;
red:=[];blue:=[];
for i from 1 to n do
for j from 1 to n do
if AB[i,j]=0 then
red:=[op(red), [i,j]];
else
blue:=[op(blue),[i,j]];
end if;
end do;
end do;
pred:=pointplot(red,symbol=circle,symbolsize=20,color=RED,axes=BOXED):
pblue:=pointplot(blue,symbol=circle,symbolsize=20,color=BLUE,axes=BOXED):
display(pred,pblue);
end proc:

> AB_display(AB);

> ij_energy:=proc(i,j)
global omega,AB,n;
local im,ip,jm,jp,ij_energy;
im:=i-1;ip:=i+1;
jm:=j-1;jp:=j+1;
if im=0 then im:=n; end if;
if ip=n+1 then ip:=1 end if;
if jm=0 then jm:=n; end if;
if jp=n+1 then jp:=1 end if;
ij_energy:=0;
if AB[i][j]<>AB[im][j] then ij_energy:=ij_energy+omega; end if;
if AB[i][j]<>AB[ip][j] then ij_energy:=ij_energy+omega; end if;
if AB[i][j]<>AB[i][jm] then ij_energy:=ij_energy+omega; end if;
if AB[i][j]<>AB[i][jp] then ij_energy:=ij_energy+omega; end if;
eval(ij_energy);
end proc:

> total_E:=0;
for i from 1 to n do
for j from 1 to n do
total_E:=total_E+ij_energy(i,j);
end do; end do;
total_E;

> sel_ij:=rand(1..n):
Energy:=[];
AB_move:=[AB_display(AB)]:
Temp:=2.0;

> for iter from 1 to 3000 do
i1:=sel_ij();
j1:=sel_ij();
i2:=sel_ij();
j2:=sel_ij();
if (AB[i1,j1]=AB[i2,j2]) then next; end if;
del_E:=-ij_energy(i1,j1)-ij_energy(i2,j2);
tmp:=AB[i1,j1];
AB[i1,j1]:=AB[i2,j2];
AB[i2,j2]:=tmp;
del_E:=del_E+ij_energy(i1,j1)+ij_energy(i2,j2);
#if del_E>0 then
if evalf(exp(-del_E/Temp))<evalf(rand()/10^12) then
tmp:=AB[i1][j1];
AB[i1,j1]:=AB[i2,j2];
AB[i2,j2]:=tmp;
else
total_E:=total_E+del_E;
Energy:=[op(Energy),total_E];
# AB_move:=[op(AB_move),AB_display(AB)]:
end if;
end do:

> total_E;

> #AB_display(AB);

> listplot(Energy);

> display(AB_move[1]);

> AB_display(AB);