1
$\begingroup$

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$$

Any advice is welcome.

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

enter image description here

Smoothed Heat Capacity

Expected behavior is like that of $x\times$Gaussian.

$\endgroup$
  • $\begingroup$ Can you add a plot of the energy vs time for a (say, intermediate) temperature? Is it possible you are not simulating long enough? $\endgroup$ – lr1985 Dec 1 '19 at 16:51
  • $\begingroup$ 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. $\endgroup$ – MurderOfCrows Dec 1 '19 at 18:25
  • $\begingroup$ @lr1985 I forgot to tag you in the reply. $\endgroup$ – MurderOfCrows Dec 1 '19 at 19:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.