Problem 2: Gibbs Sampler
Use only Python standard library modules, numpy, and matplotlib for this problem.
A Gibbs sampler generates samples from a joint probability distribution by iteratively sampling each variable conditioned on the others. In this problem, you will implement Gibbs sampling to analyze the 2-dimensional Ising model, which represents a lattice of binary spins interacting with nearest neighbors. This model is used in statistical physics to study phase transitions and collective behavior arising from local spin interactions.
Implement a Gibbs sampler for a \(N = 20\) grid (i.e. \(20 \times 20\)) where each site has state \(s_i \in \{-1, +1\}\). The system energy is: \[E(\mathbf{s}) = -\sum_{\langle i,j \rangle} s_i s_j\] where the sum is over nearest-neighbor pairs. The probability of a configuration is: \[P(\mathbf{s}) \propto \exp(-\beta E(\mathbf{s}))\] Initialize the grid randomly and store the initial state.
Run the sampler for \(\beta = \{0.2, 0.4, 0.6\}\). For each site update: \[P(\text{flip}) = \frac{1}{1 + \exp(\beta\Delta E_i)},\] where \(\Delta E_i = \text{energy change if site } i \text{ flips}\). Record the energy after every 10 iterations for 1000 total iterations. Create a figure showing \(E(\mathbf{s})\) versus iteration number for each \(\beta\).
Calculate the average magnetization: \[M = \frac{1}{N^2} \sum_{\langle i,j \rangle} s_i.\] Plot \(M\) versus iteration number. Create visualizations of the grid state at iterations \(\{0, 100, 500, 1000\}\) using:
plt.imshow(grid, cmap='binary') plt.colorbar() plt.savefig(f'state_{beta}_{iter}.png')