1D-Ising-Model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1D Ising Model\n",
    "\n",
    "This notebook will guide you through some aspects of the one-dimensional Ising model. The Ising Model originated from magnetism, but has many other applications as well. It models the interactions between neighbooring sites as\n",
    "\n",
    "\n",
    "$$\n",
    "H(\\{\\sigma\\}) = \\sum_{i,j} J \\cdot \\sigma_i\\sigma_j - \\mu \\sigma_i \\quad \\sigma_i\\in\\{-1,+1\\}\n",
    "$$\n",
    "\n",
    "Each site can be in one of two states (e.g., spin up $\\uparrow$ and spin down $\\downarrow$). If the exchange energy $J$ is positive, anti-aligned spin pairs ($\\uparrow\\downarrow$) are favoured, otherwise aligned spin pairs ($\\uparrow\\uparrow$) are energetically advantageous. The point term represents interactions with an external field $\\mu$. The spin variables $\\sigma_i$ can take on two values: $+1$ and $-1$, respectively. These choices to designate the spin states are not unique, but they make some of the math simpler."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Helping functions\n",
    "\n",
    "Before turning to solution schemes for the Ising model, we load needed libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import time\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and define some helpful functions.\n",
    "\n",
    "### Lattice state iterator\n",
    "\n",
    "For some of the solution approaches, we need to iterate over all possible states of the lattice. This is best done with an iterator. Itertors are a handy feature of python that allow to generate sequences of values that are returned by `yield` each time the iterator is accessed.\n",
    "\n",
    "This iterator returns an array containing spin variables $[\\pm 1,\\pm 1,\\dots,\\pm 1]$. \n",
    "\n",
    "> The implementation utilised the bit representation of integers, e.g., 3 -> 00000010, which allows to iterate through all possible states by simply incrementing a number.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Create a 1D lattice of spin variables\n",
    "def lattice_iterator(N):\n",
    "    '''\n",
    "    Create an iterator that iterates over all possible states {-1,+1} of a set of spin variables\n",
    "    \n",
    "    This simple implementation expoits binary representations of intergers with some padding\n",
    "    '''\n",
    "    \n",
    "    n=-1\n",
    "    M=2**N\n",
    "    d = np.zeros(N)-1;\n",
    "    \n",
    "    while n<(M-1):\n",
    "        n=n+1\n",
    "        bit = [1 if digit == '1' else -1 for digit in bin(n)[2:]]\n",
    "        d[N-len(bit):] = bit\n",
    "        yield d, n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Operators\n",
    "\n",
    "The matrix transform formalism explained below utilises **operators** acting on a function space. In computational studies, these operators take the form of matrices (similar to the formalism that you have already seen in the FEM method). Generally, an operator matrix contains the value of the operator for two states that are combined with each other through the operator in each of its matrix cells.\n",
    "\n",
    "### Magnetisation\n",
    "\n",
    "The Magnetisation operator only depends the state at hand. Drawing another analogy to the FEM formalism, it is similar to a linear form (while the Hamiltonian that we look at next is similar to a bilinear form). The matrix is of size $2^M\\times 2^M$, because this is the numer of possible states for a chain of length $M$. The magnetisation operator holds the total spin on its main diagonal for each possible spin configuration, which is constructed using the `lattice_iterator` that has been defined earlier. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def magnetisation(M):\n",
    "    m=np.zeros((2**M,2**M))\n",
    "    for s,i in lattice_iterator(M):\n",
    "        m[i,i] = np.sum(s)\n",
    "    return np.matrix(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hamiltonian\n",
    "\n",
    "The Hamiltonian is a sum over point and pair terms. The `point` function sums all spin values over an array of spins. The `pair` function sums products of spin variables $\\sigma_a\\sigma_b$ over pairs of neighbouring lattices. The benefit of this approach will become apparent later. It can be thought of as \"summing over the rugs of a ladder\".\n",
    "\n",
    "The `hamiltonian` itself computes the energy of all possible arrangements of two spin chains, and returns the resulting energies as a matrix. For, example the Hamiltonian matrix with just one spin per chain looks like this:\n",
    "\n",
    "$$\n",
    "\\begin{array}{r|cc}\n",
    "           & \\uparrow             & \\downarrow\\\\\n",
    "  \\hline\n",
    "  \\uparrow & J(+1)(+1) - \\mu(+1) & J(+1)(-1) - \\mu(-1)\\\\\n",
    "  \\downarrow & J(+1)(-1) - \\mu(+1) & J(-1)(-1) - \\mu(-1)\\\\\n",
    "\\end{array}\n",
    "\\quad\\rightarrow\\quad{\\bf H}_1 = \\left[\n",
    "\\begin{array}{cc}\n",
    "J-\\mu & -J+\\mu \\\\\n",
    "-J-\\mu & J+\\mu\\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "In this representation, spin states are described by a state vector $|\\sigma\\rangle$ (using Dirac notation, which we will learn more about later). For now, you can think of it as a fancy way of writing a vector:\n",
    "\n",
    "$$\n",
    "\\begin{array}{rcl}\n",
    "|\\uparrow\\rangle & = & [1,0]\\\\\n",
    "|\\downarrow\\rangle & = & [0,1]\\\\\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "The elements of the Hamiltonian, are then given by $H_{ij} = \\langle\\sigma_i|H|\\sigma_j\\rangle$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def point(lat,na):\n",
    "    '''\n",
    "    Sum the spins over the lattice\n",
    "    '''\n",
    "    res = 0\n",
    "    for i in range(na):\n",
    "        res += lat[i]\n",
    "\n",
    "    return res\n",
    "\n",
    "## Define functions for summing over lattice states for \n",
    "## intra-layer and trans-layer interactions\n",
    "def pair(lat_a,lat_b,n):\n",
    "    '''\n",
    "    Sum intra-layer interactions over a lattice state\n",
    "    '''\n",
    "    res = 0\n",
    "    for i in range(n):\n",
    "        res += lat_a[i]*lat_b[i]\n",
    "    return res\n",
    "\n",
    "def hamiltonian(J,mu,M):\n",
    "    H=np.zeros((2**M,2**M))\n",
    "    for s,i in lattice_iterator(M):\n",
    "        for t,j in lattice_iterator(M):\n",
    "            H[i,j] = J*pair(s,t,M) - mu*point(s,M)\n",
    "    return H"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions\n",
    "\n",
    "### Brute force\n",
    "\n",
    "The aim is to compute the thermoodynamic properties of the spin chain, namely the free energy and average spin as a function of temperature and the model parameters. Statistical Mechanics provides the tools to establish this link.\n",
    "\n",
    "The central quantity in Statistical Mechanics is the **Partition Function**:\n",
    "\n",
    "$$\n",
    "Z = \\sum_{\\{\\sigma\\}}\\exp(-\\beta\\cdot H(\\{\\sigma\\}))\n",
    "$$\n",
    "\n",
    "The partition function sums over *all possible* states $\\{\\sigma\\}$ of the spin chain, weighting each state by the exponential of the states energy $H(\\{\\sigma\\})$ divided by $kT$. It is common to use $\\beta = 1/kT$ to save some writing. This quickly becomes computationally *very* expensive. The 1D spin-chain with $N$ sites has $2^N$ possible configurations. To make this explicit, a chain with 4 sites has 16 states. A chain with 10 sites has 1024 sites, and a chain with 20 sites has 1,048,576 states.\n",
    "\n",
    "The free energy $g$ is directly computed from the partition function using the thermodynamic identity\n",
    "\n",
    "$$\n",
    "g = -kT\\cdot\\ln(Z)/N\n",
    "$$\n",
    "\n",
    "The internal energy and magnitisation is computed using explicit summation:\n",
    "\n",
    "$$\n",
    "u = \\langle H\\rangle = \\sum_{\\{\\sigma\\}}H_i\\cdot\\exp\\{-\\beta\\cdot H_i\\}\\\\\n",
    "m = \\langle \\sigma\\rangle = \\sum_{\\{\\sigma\\}}\\left(\\sum\\sigma\\right)\\cdot\\exp\\{-\\beta\\cdot H_i\\}\n",
    "$$\n",
    "\n",
    "where $\\left(\\sum\\sigma\\right)$ is the sum of spins for a given configuration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def brute_force(M,J,mu,kTs):\n",
    "    '''\n",
    "    Solve the one-dimensional Ising model by explicitly summing over all states to compute\n",
    "    partition function, free energy, internal energy, and magnetisation.\n",
    "    \n",
    "    The function returns the results as a Pandas dataframe.\n",
    "    \n",
    "    Parameters:\n",
    "        \n",
    "        M   - length of spin chain\n",
    "        J   - exchange interaction\n",
    "        mu  - external field (same units as J)\n",
    "        kTs - array with temperatures (in units of J)\n",
    "    \n",
    "    Returns:\n",
    "    \n",
    "        A Pandas dataframe with free energy g, internal energy u, magnetisation m\n",
    "        \n",
    "    Usage: \n",
    "        \n",
    "        brute_force(10,-1,0.1,np.linspace(0.1,4,30))\n",
    "    '''\n",
    "    \n",
    "    # Create a new data-frame to take results\n",
    "    df = pd.DataFrame(columns=['kT','g','u','m'],dtype=object)\n",
    "    \n",
    "    # Iterate over temperatures\n",
    "    for i,kT in enumerate(kTs):\n",
    "        beta =1/kT\n",
    "        Z = m = g = u = 0\n",
    "        \n",
    "        # Iterate over all possible spin states\n",
    "        # of a spin chain with length M\n",
    "        for l,j in lattice_iterator(M):\n",
    "            H = J*pair(l,np.roll(l,1),M)-mu*point(l,M)\n",
    "            Z += np.exp(-beta*H)\n",
    "            m += np.sum(l)*np.exp(-beta*H)\n",
    "            u += H*np.exp(-beta*H)\n",
    "\n",
    "        # Normalise sums and compute free energy\n",
    "        m = m/Z/M\n",
    "        u = u/Z/M\n",
    "        g = -1/beta*np.log(Z)/M\n",
    "        \n",
    "        # Append result to dataframe\n",
    "        s = pd.Series({'g': g, 'u': u, 'kT': kT, 'm': m, 'mu': mu, 'J': J, 'M':M}, name=i)\n",
    "        df = df.append(s)\n",
    "        \n",
    "    \n",
    "    \n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Execise 1: AFM solution of the 1D Ising model\n",
    "\n",
    "Let's investigate a simple 1D Ising chain with direct summation of the partition function. We will look at the anti-ferro-magnetic solution by setting $J<0$ and vary the "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "M = 6      # Size of the spin chain \n",
    "J = 1      # Exchange interaction\n",
    "kT = 0.1   # Temperature\n",
    "\n",
    "# Prepare plots\n",
    "df = pd.DataFrame()\n",
    "\n",
    "# Brute-Force (green)\n",
    "for mu in np.linspace(-5,5,100):\n",
    "    bf = brute_force(M,J,mu,[kT])\n",
    "    df = df.append(bf)\n",
    "\n",
    "plt.figure(figsize=(10,10));\n",
    "plt.plot(df.mu,df.m);\n",
    "plt.xlabel('External field strength $\\mu$ [eV]');\n",
    "plt.ylabel('Average magnetisation $m$ [-]');\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transfer Matrix\n",
    "\n",
    "The Transfer Matrix approach exploits a *duality* between the 1D Ising Model and a 0D quantum mechanical problem. The problem is recast into an eigenvalue problem of a $2\\times 2$ matrix, which is very efficient to solve (even analytically as we will see below). \n",
    "\n",
    "Lets start with the usual partition function\n",
    "\n",
    "$$\n",
    "Z= \\sum_{\\{\\sigma\\}} \\exp(-\\beta H_i) = \\sum_{\\{\\sigma\\}} \\exp\\left(-\\beta \\sum_j(J\\sigma_j\\sigma_{j+1} +\\mu\\sigma_j) \\right)=\\sum_{\\{\\sigma\\}}\\prod_j\\exp(-\\beta(\\underbrace{J\\sigma_j\\sigma_{j+1} +\\mu\\sigma_j}_{=H_{j,j+1}}))\n",
    "$$\n",
    "\n",
    "Now we switch to the matrix representation that is described earlier by letting\n",
    "\n",
    "$$\n",
    "T_{j,j+1} = \\exp(-\\beta\\cdot H_{j,j+1}) = \\rightarrow \\langle j|T| j+1\\rangle = \\langle j|\\exp(-\\beta H)| j+1\\rangle\n",
    "$$\n",
    "\n",
    "With this definition, the partition function becomes\n",
    "\n",
    "$$\n",
    "Z = \\sum_{\\{\\sigma\\}}\\langle 1|T| 2\\rangle\\langle 2|T| 3\\rangle\\dots\\langle N|T|1\\rangle\n",
    "$$\n",
    "\n",
    "where the last term exploits periodic boundary conditions. Now we recall that $|i\\rangle$ really only is a vector and $\\langle i|$ is its (conjugate) transpose. So $|i\\rangle\\langle i|$ is the *outer* product (or tensor product) of a unit vector. In quantum mechanics, it is called the projection operator. An important property of the projection operator is that \n",
    "\n",
    "$$\n",
    "I = \\sum_n |n\\rangle\\langle n|\n",
    "$$\n",
    "\n",
    "the sum over all projections is the identity. This can be understood as taking a vector and projecting out all its components and then summing up the components, which restores the vector. \n",
    "\n",
    "Because the sum of the partition function is over all spin configurations, it contains sums over all projection operators apart from $|1\\rangle\\langle 1|$. This greatly simplifies the sum over products, as can be demonstrated with a simple two site chain:\n",
    "\n",
    "$$\n",
    "\\begin{array}{rcll}\n",
    "Z & = & \\langle1_+|T|2_+\\rangle\\langle2_+|T|1_+\\rangle & (+,+)\\\\\n",
    "  & + & \\langle1_+|T|2_-\\rangle\\langle2_-|T|1_+\\rangle & (+,-)\\\\\n",
    "  & + & \\langle1_-|T|2_+\\rangle\\langle2_+|T|1_-\\rangle & (-,+)\\\\\n",
    "  & + & \\langle1_-|T|2_-\\rangle\\langle2_-|T|1_-\\rangle & (-,-)\\\\\n",
    "  & = & \\langle1_+|T|\\left[\\underbrace{|2_+\\rangle\\langle2_+|+|2_-\\rangle\\langle2_-|}_{=I}\\right]|T|1_+\\rangle\\\\\n",
    "  & + & \\langle1_-|T|\\left[\\underbrace{|2_+\\rangle\\langle2_+|+|2_-\\rangle\\langle2_-|}_{=I}\\right]|T|1_-\\rangle\\\\\n",
    "  & = & \\langle1_+|T^2|1_+\\rangle + \\langle1_-|T^2|1_-\\rangle\\\\\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "In general you have\n",
    "\n",
    "$$\n",
    "\\begin{array}{rcl}\\displaystyle\n",
    "Z & = & \\sum_{\\{\\sigma\\}}\\langle 1|T| 2\\rangle\\langle 2|T| 3\\rangle\\dots\\langle N|T|1\\rangle\\\\\n",
    "  & = & \\sum_{\\sigma = \\pm 1}\\langle1|T^N|1\\rangle\\\\\n",
    "  & = & {\\rm Tr}({\\bf T^N})\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "The partition function can, therefore, be computed from the trace of the N-th power of the transfer matrix through this interesting duality with quantum mechanics.\n",
    "\n",
    "Taking the $N$-th power of a matrix, however, is still an expensive operation if $N$ is large. This can be avoided by rotating the basis states so that ${\\bf T}$ becomes diagonal. In this basis, the eigenvalues $\\lambda_i$ of ${\\bf T}$ occupy the main diagonal and the trace becomes:\n",
    "\n",
    "$$\n",
    "{\\rm Tr}({\\bf T}^N) = \\sum \\lambda_i^N\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def transfer_matrix(M,J,mu,kTs):\n",
    "    '''\n",
    "    Solve the one-dimensional Ising model by computing the partition function through the \n",
    "    transfer matrix formalism.\n",
    "    \n",
    "    The function returns the results as a Pandas dataframe containing the free energy g,\n",
    "    the internal energy u, and the magnetisation m\n",
    "    \n",
    "    Parameters:\n",
    "        \n",
    "        M   - length of spin chain\n",
    "        J   - exchange interaction\n",
    "        mu  - external field (same units as J)\n",
    "        kTs - array with temperatures (in units of J)\n",
    "    \n",
    "    Returns:\n",
    "    \n",
    "        A Pandas dataframe with free energy g, internal energy u, magnetisation m\n",
    "        \n",
    "    Usage: \n",
    "        \n",
    "        transfer_matrix(10,-1,0.1,np.linspace(0.1,4,30))\n",
    "    '''\n",
    "\n",
    "    # Create a new data-frame to take results\n",
    "    df = pd.DataFrame(columns=['kT','g','u','m'],dtype=object)\n",
    "    \n",
    "    # Iterate over all tempertures\n",
    "    for i,kT in enumerate(kTs):\n",
    "        # Initialise varialbes\n",
    "        Z = m = u = 0\n",
    "        beta = 1/kT\n",
    "\n",
    "        # Compute the Hamiltonian and Transfer Matrix\n",
    "        H = hamiltonian(J,mu,1)\n",
    "        T = np.exp(-beta*H)\n",
    "        \n",
    "        # Compute eigenvalues of T and compute its eigenvalues\n",
    "        v,l = np.linalg.eig(T)\n",
    "        L = np.matrix(np.diag(v)**M)\n",
    "        Z = np.sum(v**M)\n",
    "\n",
    "        # Rotate the the Pauli matrix into the eigenbasis of T\n",
    "        U  = np.matrix(l)\n",
    "        Ui = np.matrix(np.linalg.inv(U))\n",
    "        Mg = Ui*magnetisation(1)*U\n",
    "                \n",
    "        # Compute thermodynamic quantities\n",
    "        m = np.trace(L*Mg)/Z\n",
    "        \n",
    "        g = -1/beta*np.log(Z)/M\n",
    "        \n",
    "         # Append result to dataframe\n",
    "        s = pd.Series({'g': g, 'kT': kT, 'm': m, 'mu': mu, 'J':J, 'M':m}, name=i)\n",
    "        df = df.append(s)\n",
    "    \n",
    "    # The internal energy is not easily obtained in the transfer matrix formalism,\n",
    "    # we therefore use U = -dZ/dbeta to calculate the internal energy numerically\n",
    "    if len(kTs) > 1:\n",
    "        lZ = (1/kTs)*df.g\n",
    "        df['u'] = np.gradient(lZ,1/kTs)\n",
    "    \n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2: Magnetisation of odd AFM chains\n",
    "\n",
    "The Transfer Matrix method is computationally more efficient, which allows to study longer chains. The below cell calculates the magnetisation for a longer chain with an odd number of spins."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "M = 19      # Size of the spin chain \n",
    "J = 1      # Exchange interaction\n",
    "kT = 0.1   # Temperature\n",
    "\n",
    "# Prepare plots\n",
    "df = pd.DataFrame()\n",
    "\n",
    "# Brute-Force (green)\n",
    "for mu in np.linspace(-5,5,100):\n",
    "    bf = transfer_matrix(M,J,mu,[kT])\n",
    "    df = df.append(bf)\n",
    "\n",
    "plt.figure(figsize=(10,10));\n",
    "plt.plot(df.mu,df.m);\n",
    "plt.xlabel('External field strength $\\mu$ [eV]');\n",
    "plt.ylabel('Average magnetisation $m$ [-]');\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analytical solution\n",
    "\n",
    "The transfer matrix can be solved analytically. The eigenvalues of the transfer matrix\n",
    "\n",
    "$$\n",
    "{\\bf T} = \\left[\n",
    "\\begin{array}{cc}\n",
    "\\exp(J-\\mu) & \\exp(-J+\\mu) \\\\\n",
    "\\exp(-J-\\mu) & \\exp(J+\\mu) \\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "are\n",
    "\n",
    "$$\n",
    "\\lambda_\\pm = \\exp(-\\beta J)\\cdot\\cosh(\\beta\\cdot\\mu) \\pm \\sqrt{\n",
    "    \\exp(-2\\beta J)\\cdot\\sinh^2(\\beta\\cdot\\mu) + \\exp(2\\beta/ J)\n",
    "}\n",
    "$$\n",
    "\n",
    "Once the eigenvalues are known, the partition function becomes \n",
    "\n",
    "$$\n",
    "Z = \\lambda^N_+ + \\lambda^N_- \\quad {N\\rightarrow\\infty} \\quad Z \\approx \\lambda^N_+\n",
    "$$\n",
    "\n",
    "In the thermodynamic limit ($N\\rightarrow\\infty$) the largest eigenvalue dominates and the thermodynamic properties can be computed from\n",
    "\n",
    "$$\n",
    "g = -kT\\cdot\\ln(\\lambda_+)/N\\quad\n",
    "u = -\\frac{\\partial \\ln(\\lambda_+)}{\\partial \\beta}/N\\quad\n",
    "m = \\frac{\\partial g}{\\partial \\mu}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def analytical_solution(J,mu,kTs):\n",
    "    '''\n",
    "    Solve the one-dimensional Ising model analytically through the \n",
    "    transfer matrix formalism.\n",
    "    \n",
    "    The function returns the results as a Pandas dataframe containing the free energy g,\n",
    "    the internal energy u, and the magnetisation m\n",
    "    \n",
    "    Parameters:\n",
    "        \n",
    "        J   - exchange interaction\n",
    "        mu  - external field (same units as J)\n",
    "        kTs - array with temperatures (in units of J)\n",
    "    \n",
    "    Returns:\n",
    "    \n",
    "        A Pandas dataframe with free energy g, internal energy u, magnetisation m\n",
    "        \n",
    "    Usage: \n",
    "        \n",
    "        analytical_solution(-1,0.1,np.linspace(0.1,4,30))\n",
    "    '''\n",
    "    \n",
    "    betas = 1/kTs\n",
    "\n",
    "    # Create a new data-frame to take results\n",
    "    df = pd.DataFrame(columns=['kT','g','u','m'],dtype=object)\n",
    "    \n",
    "    ## Eigenvalue\n",
    "    lmbda = np.exp(-betas*J)*np.cosh(betas*mu) \\\n",
    "          + np.sqrt(np.exp(-2*betas*J)*np.sinh(betas*mu)**2 + np.exp(2*betas*J))\n",
    "    \n",
    "    # Free Energy\n",
    "    g = -1/betas*np.log(lmbda)\n",
    "\n",
    "    # Internal energy (fairly long expression)\n",
    "    ua = np.exp(J*betas)*np.cosh(betas*mu)\n",
    "    ub = np.exp(J*betas)*mu*np.sinh(betas*mu)\n",
    "    uc = mu*np.cosh(betas*mu)*np.sinh(betas*mu) + J*(np.exp(4*J*betas)-np.sinh(betas*mu)**2)\n",
    "    ud = np.sqrt(np.exp(-2*betas*J)*(np.exp(4*J*betas)+np.sinh(betas*mu)**2))\n",
    "    ue = np.exp(-J*betas)*np.cosh(betas*mu)\n",
    "    u = -np.exp(-2*J*betas)*(-J*ua+ub+uc/ud)/(ue+ud)\n",
    "\n",
    "    # Magnetisation\n",
    "    m = np.exp(-betas*J)*np.sinh(betas*mu)/np.sqrt(np.exp(-2*betas*J)*np.sinh(betas*mu)**2 + np.exp(2*betas*J))\n",
    "\n",
    "    df['u'] = u\n",
    "    df['g'] = g\n",
    "    df['m'] = m\n",
    "    df['kT'] = kTs\n",
    "    # df['mu'] = np.ones(len(mus))*mus\n",
    "    return df\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Excercise 3: Scaling analysis\n",
    "\n",
    "Let's first check that all three methods yield similar results and then see how the \"Brute Force\" method scales with chain length compared to the \"Transfer Matrix\" method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "kTs =np.linspace(0.1,4,40)\n",
    "M = 10\n",
    "J = -1\n",
    "mu = 0.1\n",
    "\n",
    "# Prepare plots\n",
    "fig = plt.figure(figsize=(15,6))\n",
    "ax1 = plt.subplot(1,3,1)\n",
    "ax2 = plt.subplot(1,3,2)\n",
    "ax3 = plt.subplot(1,3,3)\n",
    "\n",
    "# Analytical solution (red)\n",
    "an = analytical_solution(J,mu,kTs)\n",
    "ax1.plot(np.abs(an.kT/J),an.g/np.abs(J),'r')\n",
    "ax2.plot(np.abs(an.kT/J),an.u/np.abs(J),'r')\n",
    "ax3.plot(np.abs(an.kT/J),an.m,'r')\n",
    "\n",
    "# Brute-Force (green)\n",
    "start = time.time()\n",
    "bf = brute_force(M,J,mu,kTs)\n",
    "end = time.time()    \n",
    "print(\"Brute force solution took {} sec\".format(end-start))\n",
    "ax1.plot(np.abs(bf.kT/J),bf.g/np.abs(J),'g-')\n",
    "ax2.plot(np.abs(bf.kT/J),bf.u/np.abs(J),'g-')\n",
    "ax3.plot(np.abs(bf.kT/J),bf.m,'g-')\n",
    "print(\"Avg. error Brute Force free energy       : {}\".format(np.average(np.abs(bf.g-an.g))))\n",
    "print(\"Avg. error Brute Force magnetisation     : {}\".format(np.average(np.abs(bf.m-an.m))))\n",
    "\n",
    "    \n",
    "# Transfer Matrix (blue)\n",
    "start = time.time()\n",
    "tm = transfer_matrix(M,J,mu,kTs)\n",
    "end = time.time()    \n",
    "print(\"Transfer matrix solution took {} sec\".format(end-start))\n",
    "ax1.plot(np.abs(tm.kT/J),tm.g/np.abs(J),'b-')\n",
    "ax2.plot(np.abs(tm.kT/J),tm.u/np.abs(J),'b-')\n",
    "ax3.plot(np.abs(tm.kT/J),tm.m,'b-')\n",
    "print(\"Avg. error Transfer Matrix free energy   : {}\".format(np.average(np.abs(tm.g-an.g))))\n",
    "print(\"Avg. error Transfer Matrix magnetisation : {}\".format(np.average(np.abs(tm.m-an.m))))\n",
    "\n",
    "ax1.set_xlabel('temperature kT/J [-]')\n",
    "ax1.set_ylabel('free energy g/J [-]')\n",
    "ax2.set_xlabel('temperature kT/J [-]')\n",
    "ax2.set_ylabel('internal energy u/J [-]')\n",
    "ax3.set_xlabel('temperature kT/J [-]')\n",
    "ax3.set_ylabel('magnetisation m [-]')\n",
    "\n",
    "ax2.set_xlim([0,4])\n",
    "ax3.set_ylim([-0.1,1.1])\n",
    "fig.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}