```
In [8]:
```%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
plt.style.use('bmh')
from numba import jit

```
In [23]:
```# The functions below can be used to initialise a 2D
# grid of spins, and then let it evolve towards a more
# likely configuration. See further details in
# http://web.phys.ntnu.no/~ingves/Teaching/TFY4235/Download/TFY4235_Slides_2015.pdf
# under Metropolis Monte-Carlo algorithm.
def initialise(N):
# This function initialises a grid of size N x N
# where each number is either 1 or -1
# (corresponding to spin up or down)
spins = np.ones((N, N), dtype = np.int32)
for i in range(N):
for j in range(N):
r = np.random.random()
if r < 0.5:
spins[i,j] = -1
return spins
@jit(nopython = True)
def dH(i, j, spins):
# This function calculates the change in energy
# from flipping a spin
Nx = spins.shape[0]
Ny = spins.shape[1]
J = 1
s = -spins[i,j]
return -J * s * (
spins[i-1, j] + spins[(i+1) % Nx, j] +
spins[i, j-1] + spins[i, (j+1) % Ny] )
# This function receives a considerable speedup
# from just-in-time compilation with numba, as
# it contains a for-loop. The effect of jit-compiling
# the function dH above isn't as large, but it is
# required as one cannot call a non-compiled function
# from within a compiled function.
@jit(nopython = True)
def sweep(spins, T):
# This function carries out a "sweep" over the grid,
# in the sense that it tries to flip each spin once
# on average, rejecting the flip with a certain probability
# which depends on temperature and change in energy
for i in range(spins.size):
i = np.random.randint(0, spins.shape[0])
j = np.random.randint(0, spins.shape[1])
r = np.random.random()
w = np.exp(-dH(i,j,spins)/T)
if w > r:
spins[i,j] = -spins[i,j]
return spins
def average(spins):
return np.sum(spins) / spins.size
def run(Nspins, Nsweeps, T):
# This function carries out one simulation:
# Initialise a random grid
# Carry out a number of sweeps, until the system
# can be assumed to have reached equilibrium
# Then return the average spin
# (which is proportional to the bulk magnetisation)
spins = initialise(Nspins)
for i in range(Nsweeps):
spins = sweep(spins, T)
return average(spins)

```
In [24]:
```# Set up grid size, then run a bunch of simulations
# for increasing temperature
Nspins = 40
Nsweeps = 500
Temps = np.linspace(0.1, 2, 250)
Moments = []
iout = 0
for i,T in enumerate(Temps):
# Primitive progress bar
if i > iout:
print('=', end = '')
iout += len(Temps) / 50
Moments.append(run(Nspins, Nsweeps, T))

```
```

```
In [22]:
```# Plot of average spin agains temperature
# reveals a pitchfork bifurcation.
fig = plt.figure(figsize = (9,5))
plt.scatter(Temps, Moments, c = 'k', s = 5)
plt.xlabel('Temperature')
plt.ylabel('Average spin')
plt.tight_layout()

```
```

```
In [ ]:
```