using LaTeXStrings
using Plots
gr()
"""
Magnetization of a homogeneous Ising chain
with coupling `J` and external field `h`
at temperature β
"""
m_Chain(β,J,h) = begin
s = sinh(β*h)
c = cosh(β*h)
e = exp(-4β*J)
(s + s*c/sqrt(s^2+e))/(c + sqrt(s^2+e))
end
"""
Returns the coupling matrix for a homogeneous 1d Ising chain
of length `n` of unit coupling strength
"""
function makeJring(n)
Jring = diagm(ones(n-1),1) # 1s on the first off-diagonal
Jring[1,n] = 1 # PBC
return Jring = (Jring+Jring'); # symmetrize
end
"""
Construct coupling matrix of unit strength
for a 2d square lattice with
open boundary conditions
of size `nxn`
"""
function makeJsq(n)
N = n^2
Jsq = zeros(N,N);
for k in 1:N
if k+1<=N && k%n>0
Jsq[k,k+1] = 1.
end
if k+n<=N
Jsq[k,k+n] = 1.
end
end
Jsq = Jsq+Jsq';
end
"""
Function to perform the belief propagation.
Ends when either `MaxIter` is reached or until convergence with error `err`.
Arguments:
`H`: Initial message state; get's mutated
`β`: Inverse temperature
`J`: Coupling matrix
`h`: field vector
"""
function BP!(H,β,J,h;err=1e-8,MaxIter=1000)
iter = 0
converged = false
N = size(J,1)
while iter<MaxIter && !converged
diff = 0.
for j in 1:N
# Sum over incoming messages
for i in 1:N
sum = 0.
for k in 1:N
if j==k
continue
end
sum += H[k,i]
end
if i==j
continue
end
new_hij = h[j]/N + 1/β * atanh( tanh(β*J[i,j])*tanh(β*sum))
diff = max(diff, abs(H[i,j]-new_hij))
H[i,j] = new_hij
end
end
converged = diff <= err
iter +=1
end
if !converged
println("no convergence after $MaxIter iterations!")
end
return sum(H,1)[1,:]
end
"""
Run the algorithm on a ring for a range of temperatures and
field strength `h`.
Initializes the messages itself.
Returns the magnetization over temperature.
Magnetization is taken in the middle of the ring, but for PBC this should be inconsequential.
"""
function getMagRing(n,h,βs;kw...)
hh = fill(h,n)
M = Float64[]
J = makeJring(n)
H = zeros(n,n)
H[1,2] = 1.
H[n,1] = 1.
return map(βs) do x
tanh(x*(BP!(copy(H),x,J,hh,kw...)[div(n,2)]))
end
end
"""
Run the algorithm on a square lattice for a range of temperatures
field strength `h`.
Initializes the messages to one, except the self-messages which are
set to zero.
Return the magnetization tanh(βH_i) as measured in the center of the lattice.
"""
function getMagSquare(n,h,βs;kw...)
N = n^2
hh = fill(h,N)
J = makeJsq(n)
H = ones(N,N) - eye(Float64,N)
return begin
map(βs) do β
tanh(β*(BP!(copy(H),β,J,hh;kw...)[div(N,2)]))
end
end
end
Generate data
βs = linspace(0.35,0.45,25) # temperature range
n = 10 # linear system size
hh = [0.0,] # external field strengths
## M has as columns the magnetization at a given field strength as a function of temperature
## for BP on a ring replace `getMagSquare` with `getMagRing`
M = reduce(hcat,map(h->getMagSquare(n,h,βs,MaxIter=4000), hh));
Critical temperature: between temps where magnetization is practically zero and not
βc = βs[findfirst(x->x>1e-4,M)]/2 + βs[findlast(x->x<1e-4,M)]/2
1/βc
Visualize
plot(βs,M,marker=(:auto,3),lab=["h=$h" for h in hh], lw=0, leg=:topleft)
vline!([βc], linestyle=:dash, lab=L"$ \beta_c $", lab=L"\beta_c")
xlabel!(L"\beta")
ylabel!("m")