2
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ Is MO the best place to ask for others' computer code? $\endgroup$ Jun 19, 2010 at 17:23
  • 2
    $\begingroup$ 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)" $\endgroup$ Jun 19, 2010 at 17:26
  • 3
    $\begingroup$ 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. $\endgroup$ Jun 19, 2010 at 17:39
  • $\begingroup$ 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. $\endgroup$
    – endian
    Jun 22, 2010 at 21:56

1 Answer 1

4
$\begingroup$

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;
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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