Model Isinga

Hamiltonian układu: \begin{equation*} H = -J \sum_{i, j} S_i S_j \end{equation*} S - spin
J > 0 - stała odziaływania

T = 2.26: (temp. krytyczna)

%matplotlib inline


import matplotlib.pyplot as plt
from __future__ import division
import numpy as np
from numpy.random import rand

class M_I(): 
    ## algorytm Metropolisa
    def mc_metropolis(self, konfiguracja, N, beta):
        for i in range(N):
            for j in range(N):            
                    a = np.random.randint(0, N)
                    b = np.random.randint(0, N)
                    s =  konfiguracja[a, b]
                    nb = konfiguracja[(a+1)%N,b] + konfiguracja[a,(b+1)%N] + konfiguracja[(a-1)%N,b] + konfiguracja[a,(b-1)%N]
                    cost = 2*s*nb
                    if cost < 0:	
                        s *= -1
                    elif rand() < np.exp(-cost*beta):
                        s *= -1
                    konfiguracja[a, b] = s
        return konfiguracja
    
    def wizualizacja(self, fig, konfiguracja, i, N, n_):

        X, Y = np.meshgrid(range(N), range(N))
        sp =  fig.add_subplot(3, 3, n_ )  
        plt.setp(sp.get_yticklabels(), visible=False)
        plt.setp(sp.get_xticklabels(), visible=False)      
        plt.pcolormesh(X, Y, konfiguracja, cmap=plt.cm.GnBu);
        plt.title('Czas=%d'%i)     
    plt.show()
    
    def symulacja(self):   

        N     = 100
        temp  = 2.26
        konfiguracja = 2*np.random.randint(2, size=(N,N))-1
        fig = plt.figure(figsize=(15, 15), dpi=80);    
        self.wizualizacja(fig, konfiguracja, 0, N, 1);
        
        iteracje = 1000
        for i in range(iteracje):
            self.mc_metropolis(konfiguracja, N, 1.0/temp)
            if i == 1:        self.wizualizacja(fig, konfiguracja, i, N, 2);
            if i == 10:       self.wizualizacja(fig, konfiguracja, i, N, 3);
            if i == 100:      self.wizualizacja(fig, konfiguracja, i, N, 4);
            if i == 500:     self.wizualizacja(fig, konfiguracja, i, N, 5);
            if i == 900:    self.wizualizacja(fig, konfiguracja, i, N, 6);
                
    
                 
                    
    
    
rm = M_I()

rm.symulacja()

T = 2.1: (ferromagnetyk)

%matplotlib inline


import matplotlib.pyplot as plt
from __future__ import division
import numpy as np
from numpy.random import rand

class M_I(): 
    ## algorytm Metropolisa
    def mc_metropolis(self, konfiguracja, N, beta):
        for i in range(N):
            for j in range(N):            
                    a = np.random.randint(0, N)
                    b = np.random.randint(0, N)
                    s =  konfiguracja[a, b]
                    nb = konfiguracja[(a+1)%N,b] + konfiguracja[a,(b+1)%N] + konfiguracja[(a-1)%N,b] + konfiguracja[a,(b-1)%N]
                    cost = 2*s*nb
                    if cost < 0:	
                        s *= -1
                    elif rand() < np.exp(-cost*beta):
                        s *= -1
                    konfiguracja[a, b] = s
        return konfiguracja
    
    def wizualizacja(self, fig, konfiguracja, i, N, n_):

        X, Y = np.meshgrid(range(N), range(N))
        sp =  fig.add_subplot(3, 3, n_ )  
        plt.setp(sp.get_yticklabels(), visible=False)
        plt.setp(sp.get_xticklabels(), visible=False)      
        plt.pcolormesh(X, Y, konfiguracja, cmap=plt.cm.GnBu);
        plt.title('Czas=%d'%i)     
    plt.show()
    
    def symulacja(self):   

        N     = 100
        temp  = 2.1
        konfiguracja = 2*np.random.randint(2, size=(N,N))-1
        fig = plt.figure(figsize=(15, 15), dpi=80);    
        self.wizualizacja(fig, konfiguracja, 0, N, 1);
        
        iteracje = 1000
        for i in range(iteracje):
            self.mc_metropolis(konfiguracja, N, 1.0/temp)
            if i == 1:        self.wizualizacja(fig, konfiguracja, i, N, 2);
            if i == 10:       self.wizualizacja(fig, konfiguracja, i, N, 3);
            if i == 100:      self.wizualizacja(fig, konfiguracja, i, N, 4);
            if i == 500:     self.wizualizacja(fig, konfiguracja, i, N, 5);
            if i == 900:    self.wizualizacja(fig, konfiguracja, i, N, 6);
                
    
                 
                    
    
    
rm = M_I()

rm.symulacja()

T = 2.7 (paramagnetyk):

%matplotlib inline


import matplotlib.pyplot as plt
from __future__ import division
import numpy as np
from numpy.random import rand

class M_I(): 
    ## algorytm Metropolisa
    def mc_metropolis(self, konfiguracja, N, beta):
        for i in range(N):
            for j in range(N):            
                    a = np.random.randint(0, N)
                    b = np.random.randint(0, N)
                    s =  konfiguracja[a, b]
                    nb = konfiguracja[(a+1)%N,b] + konfiguracja[a,(b+1)%N] + konfiguracja[(a-1)%N,b] + konfiguracja[a,(b-1)%N]
                    cost = 2*s*nb
                    if cost < 0:	
                        s *= -1
                    elif rand() < np.exp(-cost*beta):
                        s *= -1
                    konfiguracja[a, b] = s
        return konfiguracja
    
    def wizualizacja(self, fig, konfiguracja, i, N, n_):

        X, Y = np.meshgrid(range(N), range(N))
        sp =  fig.add_subplot(3, 3, n_ )  
        plt.setp(sp.get_yticklabels(), visible=False)
        plt.setp(sp.get_xticklabels(), visible=False)      
        plt.pcolormesh(X, Y, konfiguracja, cmap=plt.cm.GnBu);
        plt.title('Czas=%d'%i)     
    plt.show()
    
    def symulacja(self):   

        N     = 100
        temp  = 2.7
        konfiguracja = 2*np.random.randint(2, size=(N,N))-1
        fig = plt.figure(figsize=(15, 15), dpi=80);    
        self.wizualizacja(fig, konfiguracja, 0, N, 1);
        
        iteracje = 1000
        for i in range(iteracje):
            self.mc_metropolis(konfiguracja, N, 1.0/temp)
            if i == 1:        self.wizualizacja(fig, konfiguracja, i, N, 2);
            if i == 10:       self.wizualizacja(fig, konfiguracja, i, N, 3);
            if i == 100:      self.wizualizacja(fig, konfiguracja, i, N, 4);
            if i == 500:     self.wizualizacja(fig, konfiguracja, i, N, 5);
            if i == 900:    self.wizualizacja(fig, konfiguracja, i, N, 6);
                
    
                 
                    
    
    
rm = M_I()

rm.symulacja()