ising model ploting problem

J = 1;

numSpinsPerDim = 2^3;

probSpinUp = 0.5;

spin = sign(probSpinUp - rand(numSpinsPerDim, numSpinsPerDim));

kT = 1;

% Metropolis algorithm

numIters = 2^7 * numel(spin);

for iter = 1 : numIters

% Pick a random spin

linearIndex = randi(numel(spin));

[row, col] = ind2sub(size(spin), linearIndex);

% Find its nearest neighbors

above = mod(row - 1 - 1, size(spin,1)) + 1;

below = mod(row + 1 - 1, size(spin,1)) + 1;

left = mod(col - 1 - 1, size(spin,2)) + 1;

right = mod(col + 1 - 1, size(spin,2)) + 1;

neighbors = [ spin(above,col);

spin(row,left); spin(row,right);

spin(below,col)];

% Calculate energy change if this spin is flipped

dE = 2 * J * spin(row, col) * sum(neighbors);

% Boltzmann probability of flipping

prob = exp(-dE / kT);

% Spin flip condition

if dE <= 0 || rand() <= prob

spin(row, col) = - spin(row, col);

end

end

% The mean energy

sumOfNeighbors = ...

circshift(spin, [ 0 1]) ...

+ circshift(spin, [ 0 -1]) ...

+ circshift(spin, [ 1 0]) ...

+ circshift(spin, [-1 0]);

Em = - J * spin .* sumOfNeighbors;

E = 0.5 * sum(Em(:));

Emean = E / numel(spin);

% The mean magnetization

Mmean = mean(spin(:));

numSpinsPerDim = 2^3;

probSpinUp = 0.5;

J = 1;

kT = 1;

function spin = initSpins(numSpinsPerDim, probSpinUp);

spin = metropolis(spin, kT, J);

Emean = energyIsing(spin, J);

Mmean = magnetizationIsing(spin);

numSpinsPerDim = 2^3;

probSpinUp = 0.5;

J = 1;

% Temperatures to sample

numTemps = 2^9;

kTc = 2*J / log(1+sqrt(2)); % Curie temperature

kT = linspace(0, 2*kTc, numTemps);

% Preallocate to store results

Emean = zeros(size(kT));

Mmean = zeros(size(kT));

% Replace 'for' with 'parfor' to run in parallel with Parallel Computing Toolbox.

for tempIndex = 1 : numTemps

spin = initSpins(numSpinsPerDim, probSpinUp);

spin = metropolis(spin, kT(tempIndex), J);

Emean(tempIndex) = energyIsing(spin, J);

Mmean(tempIndex) = magnetizationIsing(spin);

end

figure;

plot(kT/kTc, Emean, '.');

hold on;

window = (2^-3)*numTemps - 1;

plot(kT / kTc, movmean( Emean, window));

plot(kT / kTc, movmedian(Emean, window));

hold off;

title('Mean Energy Per Spin vs Temperature');

xlabel('kT / kTc');

ylabel('<E>');

grid on;

legend('raw',...

[num2str(window) ' pt. moving mean'],...

[num2str(window) ' pt. moving median'],...

'Location', 'NorthWest');

plot(kT / kTc, Mmean, '.');

grid on;

title('Magnetization vs Temperature');

xlabel('kT / kTc');

ylabel('M');

plot(kT / kTc, abs(Mmean), '.');

hold on;

window = (2^-3)*numTemps - 1;

plot(kT / kTc, movmean( abs(Mmean), window));

plot(kT / kTc, movmedian(abs(Mmean), window));

hold off;

title('Magnetization vs Temperature');

xlabel('kT / kTc');

ylabel('|M|');

grid on;

legend('raw',...

[num2str(window) ' pt. moving mean'],...

[num2str(window) ' pt. moving median'],...

'Location', 'NorthEast');

end

why ı dont get any plots?

Answer 1

[Alan Stevens]

You have your plot commands within function "initSpins", which is never called by the coding above it!

Also check line 49, where you use "spins" while defining it!

You seem to define "numSpinsPerDim" and "probSpinUp" several times.

You haven't included functions "metropolis", "energyIsing", "magnetizationIsing", so we can't run the code to check anything even if we correct the points noted above.