Mathematica/Matlab/other for calculating Onsager's exact solution to the 2d Ising model - MathOverflow most recent 30 from http://mathoverflow.net 2013-05-19T20:09:09Z http://mathoverflow.net/feeds/question/28766 http://www.creativecommons.org/licenses/by-nc/2.5/rdf http://mathoverflow.net/questions/28766/mathematica-matlab-other-for-calculating-onsagers-exact-solution-to-the-2d-ising Mathematica/Matlab/other for calculating Onsager's exact solution to the 2d Ising model endian 2010-06-19T17:10:47Z 2010-06-19T17:30:00Z <p>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. </p> http://mathoverflow.net/questions/28766/mathematica-matlab-other-for-calculating-onsagers-exact-solution-to-the-2d-ising/28768#28768 Answer by Steve Huntsman for Mathematica/Matlab/other for calculating Onsager's exact solution to the 2d Ising model Steve Huntsman 2010-06-19T17:30:00Z 2010-06-19T17:30:00Z <p>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...</p> <hr> <pre><code>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; </code></pre>