MathOverflow is a question and answer site for professional mathematicians. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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.

share|cite|improve this question
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:…) 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
up vote 4 down vote accepted

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

% 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;

% 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);

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

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);

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

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

y = Q;
share|cite|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.