Phase transition in Ising model with ext. field

# Python3
# see: < https://physics.stackexchange.com/q/348917 >

import numpy as np
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
#==============================================================================80

#-- mean-field equation: m = tanh( Tc/T * (h + m) )

# m  := magnetic order parameter
# h  := rescaled external magnetic field = h_ext / Tc
# t  := T/Tc , rescaled temperature
# HM := Tc/T * (h + m) = (h + m)/t

#-- MF solver
def solveMF(hs, t):
    # solves MF eq for the given h-values at a rescaled temperature t

    MFeq = lambda m, h: m - np.tanh( (m+h)/t )
    ms = map(lambda h: fsolve(lambda m: MFeq(m, h), x0= -2 if h < 0 else +2), hs)

    return tuple(ms)
#-------------------

#-- thermodynamic quantities
MH = lambda m, h, t: (m + h)/t

# free energy, f(h, T; m) = F/N
FE  = lambda m, h, t: 0.5 * m**2 - t * np.log( 2*np.cosh(MH(m,h,t)) )

# dm/dt
dM_dT = lambda m, h, t: MH(m, h, t) / (1 - t/(1 - m**2))

# entropy, s(h, T; m) = S/N
def S(m, h, t):
    mh0 = MH(m, h, t)
    return np.log( 2*np.cosh(mh0) ) - m * mh0

# specific heat at const. volume, c_V(h, T; m) = C_V/N
CV = lambda m, h, t: -t * MH(m, h, t)**2 / (1 - t/(1 - m**2))

# magnetic susceptibility
Xh = lambda m, h, t: -1 / (1 - t/(1 - m**2))

#==============================================================================80

#-- free energy and entropy difference for the case where a finite opposite h
# is applied to a system prepared in an ordered phase (T < Tc)

t0 = 0.5  # T < Tc
hs = (1e-3, -1)  # initial and final magnetic fields
ms = solveMF(hs, t0)
fs = tuple( map(lambda z: FE(z[0], z[1], t0), zip(ms, hs)) )
ss = tuple( map(lambda z: S(z[0], z[1], t0) , zip(ms, hs)) )
df = fs[1]-fs[0]
ds = ss[1]-ss[0]
print("free-energy change: f(h=%.3f) - f(h=%.3f) = %.3f" % (hs[1], hs[0], df) )
print("entropy change    : s(h=%.3f) - s(h=%.3f) = %.3f" % (hs[1], hs[0], ds) )
print("heat exchange     : dq = T ds = %.3f" % (ds/t0) )

#-- plots

#--- graphical solution of MF eq
h0 = -0.1                       # rescaled magnetic field
pTitle0  = r"$T < T_c$" if t0 < 1 else r"$T > T_c$" if t0 > 1 else r"$T = T_c$"

Mmax = 2                        # max of m
ms = np.linspace(-Mmax, Mmax, 1<<10)  # m values
RHSs = list( map(lambda m_: np.tanh( (m_+h0)/t0 ), ms) )  # RHS vals.
fs = list(map(lambda m_: FE(m_, h=h0, t=t0), ms))  # free-energy vals.

plt.plot(ms, ms, color='green', label='LHS')
plt.plot(ms, RHSs, color='orange', label='RHS')
plt.plot(ms, fs, color='crimson', linestyle='-', lw = 1, label='f')

pTitle = "Graphical solution to MF eq\n" + pTitle0 + "\n" + r"$h_0$"+ "= %.3f" % h0
plt.title(pTitle, fontsize=10)
plt.xlabel("m")
plt.grid()
plt.legend()
plt.show()

#--- thermodynamic quantities
Hmax = 1                                           # max of h
hs = np.linspace(-Hmax, Hmax, 1<<10, dtype=float)  # h-field vals.
ms = solveMF(hs, t0)                               # solve for m

# m-h plot
plt.plot(hs, ms,
         color='blue', marker='.', ms=0.2, lw = 0, label=r'$m$')

plt.title(pTitle0)
plt.xlabel("h")
plt.ylabel("m")
#plt.legend()
plt.show()

# entropy
ss = list( map(lambda z: S(z[0], z[1], t0), zip(ms, hs)) )

plt.plot(hs, ss,
         color='red', lw = 1, label='s')

plt.title(pTitle0)
plt.xlabel("h")
plt.ylabel("s")
#plt.legend()
plt.show()

# magnetic susceptibility
xhs = list( map(lambda z: Xh(z[0], z[1], t0), zip(ms, hs)) )
plt.plot(hs, xhs,
         color='orange', lw = 1, label=r'$\chi_h$')

plt.title(pTitle0)
plt.xlabel("h")
plt.ylabel(r'$\chi_h$')
#plt.legend()
plt.show()

# specific heat
cvs = list( map(lambda z: CV(z[0], z[1], t0), zip(ms, hs)) )
plt.plot(hs, cvs,
         color='green', lw = 1, label=r'$c_V$')


plt.title(pTitle0)
plt.xlabel("h")
plt.ylabel(r'$c_V$')
#plt.legend()
plt.show()
#==============================================================================80