import numpy as np
import matplotlib.pyplot as plt

Le modèle d'Ising

On veut simuler une configuration $\sigma\in \{-1,1\}^{\{1,\ldots,N\}^2}$ de la mesure de Gibbs $$ \mathbb P(\sigma)=\frac1{Z_\beta}\mathrm e^{-\beta H(\sigma)} $$ à l'aide de l'algorithme Metropolis-Hastings, avec $$ H(\sigma)=-\sum_{\substack{u\sim v}}\sigma_u \sigma_v. $$

Si on propose un sommet uniformement au hasard et on change son signe, on a $$ \rho(\sigma,\tau_w\sigma)= \frac{\mathrm e^{-\beta H(\tau_w\sigma)}}{\mathrm e^{-\beta H(\sigma)}} =\mathrm e^{-\beta D(\sigma,w)} $$ avec $$ D(\sigma,w):=H(\tau_w\sigma)-H(\sigma)= 2 \sigma_w\sum_{u\sim w} \sigma_u. $$

# si w=(i,j), alors M[i,j] représente sigma(w)

def rho(M,i,j):
    D = 2*M[i,j]*( M[i-1,j] + M[(i+1)%N,j] + M[i,(j+1)%N] + M[i,j-1] ) # % veut dire "modulo"
    return np.exp(-beta*D)

def transition(M):
    (i,j) = np.random.randint(N,size=2)
    U = np.random.rand()
    if U < rho(M,i,j):
        M[i,j]*=-1 
    return M

def IsingMH(N,beta,M0,Nstep):
    M = M0
    for n in range(Nstep):
        M = transition(M)
    return M

Condition initiale : aléatoire (haute énergie)

# Conditions initiales : aléatoire

def M0rand(N): return(2*np.random.randint(2,size=(N,N))-1)

plt.imshow(M0rand(50))
plt.title("Configuration initiale (random)")
plt.show()

N = 50
# beta = 1     # different de beta = 0.1. Transition ? 
betaCritique = np.log(1+np.sqrt(2))/2
beta = betaCritique + 0.1

M0 = M0rand(N)
Nstep = int(1e6)
IsingRand =  IsingMH(N,beta,M0,Nstep)
plt.imshow(IsingRand)
plt.title("Configuration finale (random)")
plt.show()

Condition initiale : alignée (énergie minimale)

# Condition initiale : spins alignés
def M0align(N): return(np.ones((N,N)))

plt.imshow(M0align(N))
plt.title("Configuration initiale (alignée)")
plt.show()

N = 50 
beta = betaCritique - 0.1 
M0 = M0align(N)
Nstep = int(1e6)
IsingRand =  IsingMH(N,beta,M0align(N),Nstep)
plt.imshow(IsingRand)
plt.title("Configuration finale (alignée)")
plt.show()