# Mathematica/Matlab/other for calculating Onsager's exact solution to the 2d Ising model

Would anybody be able to share a Mathematica/Matlab/other script for calculating Onsager's exact solution for the magnetisation of the 2d Ising model? I would be most grateful of one in order to test my MC simulations of the system.

-
Is MO the best place to ask for others' computer code? –  Robin Chapman Jun 19 '10 at 17:23
This isn't even about code. The expression for the magnetization is simple (see, e.g. here: en.wikipedia.org/wiki/…) and the corresponding MATLAB expression would be "m = ( 1 - ( sinh(2*betaE1) * sinh(2*betaE2) )^-2 )^(1/8)" –  Steve Huntsman Jun 19 '10 at 17:26
Robin, the people most likely to be familiar with Onsager's solution to the Ising model are mathematical physicists. Consequently, yes, this is the best place to ask for such a script. –  Tom LaGatta Jun 19 '10 at 17:39
Thanks. I had asked for a script as I wanted to be certain the the 1/sinh^4 expression was completely correct. Onsager didn't quote it that way, instead it came out as a particularly nasty looking integration, and Yang just confused me even further. I shall go with the analytic expression above though, although it obviously breaks down right at the phase transition. –  endian Jun 22 '10 at 21:56

Not that it's directly relevant, but I have code for the generator matrix of a 1D Glauber-Ising model that could probably be reworked into 2D...

function y = glauber1d(symb,n,varargin);

% produces the generator matrix etc for a 1D Glauber-Ising model of n spins
% call as either glauber1d(1,n) for a less complete symbolic result, or
% glauber1d(0,5,[a,mu,H,kT,J]) for a more complete numerical result--i.e.,
% symb is a flag indicating whether or not to use symbolic calculations
% (this requires the symbolic toolbox in order to work)

% a (Glauber's alpha) is the spin flip rate, depends on the coupling
%   between the GI system and the bath;
% mu is the magnetic moment associated with the spins;
% H is the magnetic field strength;
% kT is (well, you know);
% J is the exchange energy

if symb     % SYMBOLICS
syms a b g real;
else        % NUMERICS
args = varargin{1};
a   = args(1);
mu  = args(2);
H   = args(3);
kT  = args(4);
J   = args(5);
b = tanh(mu*H/kT);   % Glauber's beta (NOT 1/kT)
g = tanh(2*J/kT);   % Glauber's gamma
end

% produce an array with rows equal to spin configurations
temp = dec2bin(0:((2^n)-1),n);
for j = 1:2^n
for k = 1:n
s(j,k) = 2*str2num(temp(j,k))-1;
end
end

% obtain spin flip rates
for j = 1:2^n
for k = 1:n
km = mod(k-2,n)+1;
kp = mod(k,n)+1;
temp = (g/2) * (b - s(j,k)) * (s(j,km) + s(j,kp));
w(j,k) = (a/2) * (1 - b*s(j,k) + temp);
end
end

% generator matrix
if symb
Q = sym(zeros(2^n));
else
Q = zeros(2^n);
end

for j1 = 1:2^n
for j2 = 1:2^n
if sum(abs( s(j1,:) - s(j2,:) )) == 2   % single spin flip
% now find out which spin gets flipped
k0 = find( s(j1,:) - s(j2,:) );
Q(j1,j2) = w(j1,k0);
end
end
end

if symb
Q = simplify( Q - diag(sum(Q,2)) );
else
Q = Q - diag(sum(Q,2));
end

% invariant distribution p (if you want it)
if 2^n - 1 - rank(Q)
'error'
y = 0;
return;
else
p0 = null(Q')';
end
if symb, simplify(p0); end
sp0 = sum(p0);
if symb, simplify(sp0); end
p = p0 / sp0;   % invariant distribution

y = Q;

-