# Using MATLAB to simulate the Ising Model

I am using MATLAB to simulate a 1D Ising Chain. I am running into an issue where when trying to find heat capacity, my system has a tremendous amount of noise. I'll post my code and an image of the heat capacity (as well as it smoothed 1000 times).

I have tried to smooth the energies first, I've tried having the "equilibrium energy" be the average of the last 5, 10, 100, 1000 timesteps (runs in the code), and nothing has worked. The max should be around 450, but smoothing will change the values.

My heat capacity definition is $$C = \dfrac{\partial E}{\partial T}$$ But I have also tried using this definition to no avail: $$C = kT \sigma_E^2$$

## Edit 1

I changed the code so that the energy is smoothed before it is differentiated, and I smoothed it with different inputs. It somewhat resolved the problem.

%% Ising1D.m is a script designed to simluate a 1d Ising Chain at different Hfields and Temps
% The Metropolis Algorithm will be used
function [k, T, J, N, E, C, M, X] = Ising1D
%% Initializing the chain and other constants
N = 1000; % Number of spins
J = 1;
k = 1;
H = linspace(0,10,100); % Magnetic Energy
T = linspace(.1,10,100); % Temperature
flip = -1; % If I want to flip a spin, I'll multiply it by this.

%% Applying Metropolis Algorithm
time = 30*N; % Number of times I run the algorithm / mcs / Monte Carlo Steps
E = zeros(length(T), length(H), time); % A vector consisting of energy of system at different points in algorithm.
M = zeros(length(T), length(H)); % Magnetization
randSpin = randi(N,time,1); % An array of randomly chosen spins
randNum = rand(time,1); % An array of random numbers between 0 and 1.
for Tindex = 1:length(T)
for Hindex = 1:length(H)
temporary = rand(1,N); % Building my chain of spins randomly
ups = temporary >= .5; % Up spins
downs = -(temporary < .5); % Down spins
chain = ups + downs; % Random initial set of spins.
site = 1:(N-1);
E(Tindex, Hindex, 1) = -J.*sum(chain(site).*chain(site+1)) - Hindex.*(sum(chain)); % Initialize energy
for run = 2:time
if randSpin(run) == 1 % calculate change in energy
dE = 2*J*chain(randSpin(run))*chain(randSpin(run)+1) + 2*Hindex*chain(randSpin(run));
elseif randSpin(run) == N
dE = 2*J*chain(randSpin(run)-1)*chain(randSpin(run)) + 2*Hindex*chain(randSpin(run));
else
dE = 2*J*chain(randSpin(run))*(chain(randSpin(run)-1) + chain(randSpin(run)+1)) +  + 2*Hindex*chain(randSpin(run));
end
if ( dE < 0 )
E(Tindex, Hindex, run) = E(Tindex, Hindex, run-1) + dE; % Update Energy.
chain(randSpin(run)) = flip*chain(randSpin(run)); % Flip the spin for good.
else
p = exp(-dE/(k*T(Tindex)));
if (randNum(run) <= p) % Accept the change
E(Tindex, Hindex, run) = E(Tindex, Hindex, run-1) + dE;
chain(randSpin(run)) = flip*chain(randSpin(run));
else % Reject the change
E(Tindex, Hindex, run) = E(Tindex, Hindex, run-1);
end
end
end
M(Tindex, Hindex) = sum(chain); % M = N*avg(mu) = N*(ups+downs)/n = sum(chain)
end
end

%% Smooth the data
E(:, 1, end) = smooth(E(:, 1, end), 0.15, 'rloess');
M(:,1) = smooth(M(:,1), .5, 'rloess');

%% Calculate and plot C
C = diff(E(:, 1, end))./diff(T)';
C = smooth( C/(N*k), .75, 'rloess'); % Scaled heat capacity
figure('Color', 'w');
subplot(1,2,1);
plot(linspace(0,8,length(C)), C); % The random energies at high temps make this very noisy. Smooth reduces noise via moving average.
title('Heat Capacity [C] vs Temperature [T]');
xlabel('T');
ylabel('C');
legend('H = 0');

%% Calculate X
X = diff(M,1,2)./diff(H);
subplot(1,2,2);
hold on;
plot(H(1:end-1), smooth(X(1,:), 1000)); plot(H(1:end-1), smooth(X(20,:), 1000)); plot(H(1:end-1), smooth(X(40,:), 1000)); plot(H(1:end-1), smooth(X(60,:), 1000)); plot(H(1:end-1), smooth(X(80,:), 1000)); plot(H(1:end-1), smooth(X(100,:), 1000));% Smoothing reduces noise but alters the data, namely by making the values smaller.
title('\chi vs Magnetic Field Strength [H]');
xlabel('H');
ylabel('\chi');
legend('T = 0.1', 'T = 2', 'T = 4', 'T = 6', 'T = 8', 'T =10');

end


## Energy vs Temp ## Not Smoothed Heat Capacity ## Smoothed Heat Capacity Expected behavior is like that of $$x\times$$Gaussian.

• Can you add a plot of the energy vs time for a (say, intermediate) temperature? Is it possible you are not simulating long enough? – lr1985 Dec 1 '19 at 16:51
• I can do that, and I'll add the updated code as well. It's possible that I'm not simulating long enough, but the energy seems to get very close to equilibrium well before the simulation ends. – MurderOfCrows Dec 1 '19 at 18:25
• @lr1985 I forgot to tag you in the reply. – MurderOfCrows Dec 1 '19 at 19:05