.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/statphys/ising/plot_metropolis_algorithm.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_statphys_ising_plot_metropolis_algorithm.py: The Metropolis algorithm: sampling Boltzmann without the partition function ============================================================================== Metropolis, the Rosenbluths and the Tellers (1953) showed how to draw states with Boltzmann probability :math:`e^{-\beta E}/Z` without ever computing :math:`Z`. Propose a small random change; accept it with probability :math:`\min(1, e^{-\beta\Delta E})`. The rule satisfies detailed balance, :math:`p_a P(a\to b) = p_b P(b\to a)`, so the chain's stationary distribution is exactly Boltzmann's, from any starting state. A 3x3 periodic Ising lattice has 512 states, few enough to compute the exact distribution. This example checks the acceptance rule, runs :func:`~physicskit.statphys.core.monte_carlo.metropolis_sweep_ising` (via :class:`~physicskit.statphys.chapters.ising_lattice.Ising2D`) from an ordered and a random start, and compares the sampled energy histogram with the exact one, at a temperature where direct sampling of random configurations would almost never hit the important states. .. GENERATED FROM PYTHON SOURCE LINES 22-31 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from physicskit.statphys.chapters.ising_lattice import Ising2D L, J, T = 3, 1.0, 2.5 beta = 1.0 / T N = L * L .. GENERATED FROM PYTHON SOURCE LINES 32-34 Exact Boltzmann distribution -------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 34-44 .. code-block:: Python states = ((np.arange(2**N)[:, None] >> np.arange(N)) & 1) * 2 - 1 spins = states.reshape(-1, L, L) E = -J * (np.sum(spins * np.roll(spins, 1, 1), axis=(1, 2)) + np.sum(spins * np.roll(spins, 1, 2), axis=(1, 2))) levels, degeneracy = np.unique(E, return_counts=True) p_exact = degeneracy * np.exp(-beta * (levels - levels.min())) p_exact /= p_exact.sum() print("E g(E) Boltzmann P(E) uniform P(E)") for e, gd, p in zip(levels, degeneracy, p_exact): print(f"{e:5.0f} {gd:6d} {p:14.5f} {gd / 2**N:12.5f}") .. rst-class:: sphx-glr-script-out .. code-block:: none E g(E) Boltzmann P(E) uniform P(E) -18 2 0.56932 0.00391 -10 18 0.20886 0.03516 -6 48 0.11245 0.09375 -2 198 0.09365 0.38672 2 144 0.01375 0.28125 6 102 0.00197 0.19922 .. GENERATED FROM PYTHON SOURCE LINES 45-52 The acceptance rule and detailed balance -------------------------------------------- For a single spin flip in the Ising model, :math:`\Delta E` can only be :math:`-8J,-4J,0,4J,8J`. Moves downhill are always taken; uphill moves only with the Boltzmann factor. The ratio of forward to backward acceptance equals the ratio of Boltzmann weights, which is detailed balance. .. GENERATED FROM PYTHON SOURCE LINES 52-56 .. code-block:: Python for dE in (-8, -4, 0, 4, 8): a_fwd, a_back = min(1, np.exp(-beta * dE)), min(1, np.exp(beta * dE)) print(f"dE = {dE:+d}J: accept {a_fwd:.4f}, reverse {a_back:.4f}, ratio {a_fwd / a_back:.4f} = exp(-beta dE) {np.exp(-beta * dE):.4f}") .. rst-class:: sphx-glr-script-out .. code-block:: none dE = -8J: accept 1.0000, reverse 0.0408, ratio 24.5325 = exp(-beta dE) 24.5325 dE = -4J: accept 1.0000, reverse 0.2019, ratio 4.9530 = exp(-beta dE) 4.9530 dE = +0J: accept 1.0000, reverse 1.0000, ratio 1.0000 = exp(-beta dE) 1.0000 dE = +4J: accept 0.2019, reverse 1.0000, ratio 0.2019 = exp(-beta dE) 0.2019 dE = +8J: accept 0.0408, reverse 1.0000, ratio 0.0408 = exp(-beta dE) 0.0408 .. GENERATED FROM PYTHON SOURCE LINES 57-62 Two starts, one distribution -------------------------------- From the all-up state and from a random state the running average of the energy converges to the same exact value, and the histogram of visited energies matches Boltzmann's distribution. .. GENERATED FROM PYTHON SOURCE LINES 62-97 .. code-block:: Python E_exact = np.sum(levels * p_exact) traces, hist = {}, {} for label, ordered in (("ordered start", True), ("random start", False)): model = Ising2D(L=L, J=J, seed=1953 + ordered) model.reset(ordered=ordered) Es = [] for _ in range(40000): model.sweep(beta) Es.append(model.energy()) Es = np.array(Es) traces[label] = np.cumsum(Es) / np.arange(1, Es.size + 1) hist[label] = np.array([np.mean(Es[1000:] == e) for e in levels]) print(f"{label}: = {Es[1000:].mean():.4f} (exact {E_exact:.4f})") fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4)) for label, tr in traces.items(): ax1.semilogx(np.arange(1, tr.size + 1), tr, label=label) ax1.axhline(E_exact, color="k", ls="--", label="exact ") ax1.set_xlabel("Metropolis sweeps") ax1.set_ylabel("running average of E") ax1.set_title("Any start converges to the same ensemble") ax1.legend(fontsize=8) w = 1.2 ax2.bar(levels - w / 2, p_exact, width=w, color="0.7", label="exact Boltzmann") ax2.bar(levels + w / 2, hist["random start"], width=w, color="firebrick", alpha=0.8, label="Metropolis") ax2.bar(levels, degeneracy / 2**N, width=0.4, color="steelblue", label="uniform random states") ax2.set_yscale("log") ax2.set_xlabel("energy E") ax2.set_ylabel("probability") ax2.set_title(f"3x3 Ising at T = {T}") ax2.legend(fontsize=8) fig.tight_layout() plt.show() .. image-sg:: /api/gallery/statphys/ising/images/sphx_glr_plot_metropolis_algorithm_001.png :alt: Any start converges to the same ensemble, 3x3 Ising at T = 2.5 :srcset: /api/gallery/statphys/ising/images/sphx_glr_plot_metropolis_algorithm_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none ordered start: = -13.0581 (exact -13.1591) random start: = -13.1855 (exact -13.1591) .. GENERATED FROM PYTHON SOURCE LINES 98-101 Check ----- Detailed balance makes Metropolis sample exp(-beta E) from any start. .. GENERATED FROM PYTHON SOURCE LINES 101-104 .. code-block:: Python assert np.sum(degeneracy) == 2**N assert np.max(np.abs(hist["random start"] - p_exact)) < 0.01 assert abs(traces["random start"][-1] - E_exact) < 0.15 and abs(traces["ordered start"][-1] - E_exact) < 0.15 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.176 seconds) .. _sphx_glr_download_api_gallery_statphys_ising_plot_metropolis_algorithm.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_metropolis_algorithm.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_metropolis_algorithm.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_metropolis_algorithm.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_