{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[],"authorship_tag":"ABX9TyNH+62mX02jnQFv+k5n6p/i"},"kernelspec":{"name":"python3","display_name":"Python 3"}},"cells":[{"cell_type":"markdown","metadata":{"id":"HwvaNgJ578Ua"},"source":["# MCMC simulation of the Ising model\n","\n","The goal of the tutorial is to implement the Metropolis algorithm to sample configurations from the Ising model on the discrete torus in dimension $d=2$."]},{"cell_type":"code","metadata":{"id":"1SICXgNjpQYO"},"source":["import numpy as np\n","import matplotlib.pyplot as plt"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"az9yXzVw9mhh"},"source":["# 1. The case $\\beta=0$\n","\n","A configuration is represented by an array of size $N \\times N$. Write a function `ising_0(N)` which returns such an array with cells containing spins independently uniformly distributed in $\\{-1,1\\}$."]},{"cell_type":"code","metadata":{"id":"13vhUcUm8qUe"},"source":["def ising_0(N):\n","  # Your code here"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"69xFZVXe_Kt2"},"source":["The function `plot_configuration(conf)` takes in argument a configuration and plots a figure in which spins $+$ are represented by a black pixel and spins $-$ are represented by a white pixel."]},{"cell_type":"code","metadata":{"id":"tmJ_1I7l-hps"},"source":["def plot_configuration(conf):\n","  fig, ax = plt.subplots()\n","  ax.matshow(conf,cmap=\"Greys\",vmin=-1,vmax=1)\n","  ax.set_axis_off()\n","\n","plot_configuration(ising_0(1000))"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"M4kaMHO_FsIZ"},"source":["# 2. Metropolis algorithm\n","\n","Write a function `sigma(conf,i,j)` which takes in argument a configuration `conf` and two indices `i` and `j` between $0$ and $N-1$, and returns the sum of the spins over the neighbours of the vertex with coordinates $(i,j)$.\n","\n","We recall that boundary conditions are periodic, so that for example, the neighbours of the vertex $(0,0)$ are $(0,N-1)$, $(0,1)$, $(N-1,0)$ and $(1,0)$. You may therefore fruitfully use the operator `%`."]},{"cell_type":"code","metadata":{"id":"1wy3Gnw8GyfY"},"source":["def sigma(conf,i,j):\n","  # Your code here"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"8d0P19C7NpdB"},"source":["Implement the Metropolis algorithm, with the rule of your choice, to return a configuration of size $N \\times N$, at inverse temperature $\\beta$, in $n_\\mathrm{max}$ iterations.\n","\n","You may take as an initial configuration one returned by `ising_0`."]},{"cell_type":"code","metadata":{"id":"mVu8IfvUF4WJ"},"source":["# Your code here"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"LfhUsv-LPCLt"},"source":["# 3. Observation of the phase transition\n","\n","For any $\\beta \\geq 0$, we set\n","$$p_\\beta = 1-\\mathrm{e}^{-2\\beta} \\in [0,1].$$\n","The critical inverse temperature $\\beta_\\mathrm{c}$ corresponds to\n","$$p_{\\beta_\\mathrm{c}} = \\frac{\\sqrt{2}}{1+\\sqrt{2}} \\simeq 0,585786.$$\n","\n","We draw typical configurations, for $N=250$, with several values of $p_\\beta$."]},{"cell_type":"code","metadata":{"id":"OvCy33eIRQAp"},"source":["val_p = [.1, .5, .6, .8]\n","\n","for p in val_p: plot_configuration(ising_MH(250,-.5*np.log(1-p),10**6))"],"execution_count":null,"outputs":[]}]}