(*in this notebook we explore the exact solution of the Ising model in 2d*)
(*-beta f= lnZ*)

k[b_,J_]:=Sinh[2 b J]^(-2)

lnZ[b_,J_]:=Log[2]/2 + 1/(2 Pi)NIntegrate[Log[Cosh[2 b J]^2+ Sqrt[1+k[b,J]^2- 2 k[b,J] Cos[2 theta]]/k[b,J] ],{theta,0,Pi}]

energy[b_,J_]:= - J Coth[2 b J] (1 + 2/Pi (2 Tanh[2 b J]^2-1) NIntegrate[1/Sqrt[1-4k[b,J](1+k[b,J])^(-2)Sin[theta]^2],{theta,0,Pi/2}])

Needs["NumericalCalculus`"]

specificheat[b_,J_]:=-b^2ND[- J Coth[2 bp J] (1 + 2/Pi (2 Tanh[2 bp J]^2-1) NIntegrate[1/Sqrt[1-4k[bp,J](1+k[bp,J])^(-2)Sin[theta]^2],{theta,0,Pi/2}]),bp,b]

m[b_,J_]:=(1-Sinh[2 b J]^-4)^(1/8)

Plot[m[1/T,1],{T,0,5},PlotRange->{0,1}]

(*the spontaneous magnetisation at h=0 is nonzero below T_c=2.2...*)

Plot[lnZ[1/T,1],{T,0,5}]

(*note the free energy does not have any peculiar feature at the phase transition at T_c*)

Plot[energy[1/T,1],{T,0,5}]

(*the energy has a large slope near T_c, as does the entropy below*)

Plot[lnZ[1/T,1]+energy[1/T,1]/T,{T,0,5}]

Plot[specificheat[1/T,1],{T,.1,5}]

(* you can find a nice collection of results at http://quantumtheory.physik.unibas.ch/people/bruder/Semesterprojekte2007/p1/index.html *)