Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
clear;
page_output_immediately(1);
% obliczanie energii lancucha
function E = energy(chain, J, H)
E = -sum(chain)*H;
E += -J*sum(chain.*[chain(end) ; chain(1:end-1)]);
endfunction
% generator LCG
function value = LCGrand()
choice = 1;
global state;
% parametry generatora LCG
a = [69069 ; 16807 ; 1664525 ; 1103515245];
c = [1 ; 0 ; 1013904223 ; 12345];
m = [ power(2,32); power(2,32) - 1; power(2,32); power(2,32)];
state = mod(a(choice)*state+c(choice),m(choice));
value = state/(m(choice)-1);
endfunction
function [E,M] = generateChain(H,J,T,N,iterations, steps, generator)
E = 0;
M = 0;
%lancuch losowych spinow
chain = round(rand(N,1)).*2-1;
for cnt1=1:steps
for cnt2=1:iterations
% losowanie nru odwracanego spinu
i = round(generator()*(N-1)+1);
% odwrocenie pojedynczego spinu
trial_chain = chain;
trial_chain(i) *= -1;
% roznica energii pomiedzy nowym i starym lancuchem
dE = energy(trial_chain, J, H) - energy(chain, J, H);
if(dE <= 0)
% jezeli nizsza energia nowego lancucha - odwracamy
chain = trial_chain;
else
if(exp(-dE/T) >= generator())
% jezeli energia jest wyzsza
% to z p-bstwem w/w przyjmujemy nowy lancuch
chain = trial_chain;
endif
endif
endfor
E += energy(chain,J,H)/steps;
M += abs(sum(chain))/N/steps;
endfor
endfunction
%{---------------------------------------------------------------%}
imgParam = '-S1000,800';
imgExt = 'png';
% ilosc krokow Monte Carlo
MCsteps = 50;
% ilosc iteracji do termalizacji ukladu
iterations = 1500;
% ilosc spinow
N = 100;
T = 0.0001:0.1:12;
%%
% ziarno generatora LCG
global state = 1;
count = length(T);
E = zeros(count,1);
M = zeros(count,1);
disp(':> Uruchomiono generowanie wykresow');
t0 = cputime();
%%
H = 0;
J = 1;
for i=1:count
[E(i), M(i)] = generateChain(H,J,T(i),N,iterations,MCsteps,@rand);
endfor
disp(':> Zakonczono generowanie ukladu 1.');
disp([':> Czas generacji: ' num2str((t1 = cputime()) - t0)]);
figure(1);
plot(T,E)
xlabel('Temperatura');
ylabel('Energia');
title(['Energia w funkcji temperatury J = ' num2str(J) ' H = ' num2str(H) ' generator rand()']);
grid;
print(['EJ' num2str(J) 'H' num2str(H) 'rand.' imgExt], imgParam);
close(1);
figure(2);
plot(T,M)
xlabel('Temperatura');
ylabel('Srednia magnetyzacja');
title(['Srednia magnetyzacja w funkcji temperatury J = ' num2str(J) ' H = ' num2str(H) ' generator rand()']);
grid;
print(['MJ' num2str(J) 'H' num2str(H) 'rand.' imgExt], imgParam);
close(2);
H = 0;
J = -1;
for i=1:count
[E(i), M(i)] = generateChain(H,J,T(i),N,iterations,MCsteps,@rand);
endfor
disp(':> Zakonczono generowanie ukladu 2.');
disp([':> Czas generacji: ' num2str((t2 = cputime()) - t1)]);
figure(3);
plot(T,E)
xlabel('Temperatura');
ylabel('Energia');
title(['Energia w funkcji temperatury J = ' num2str(J) ' H = ' num2str(H) ' generator rand()']);
grid;
print(['EJ' num2str(J) 'H' num2str(H) 'rand.' imgExt], imgParam);
close(3);
figure(4);
plot(T,M)
xlabel('Temperatura');
ylabel('Srednia magnetyzacja');
title(['Srednia magnetyzacja w funkcji temperatury J = ' num2str(J) ' H = ' num2str(H) ' generator rand()']);
grid;
print(['MJ' num2str(J) 'H' num2str(H) 'rand.' imgExt], imgParam);
close(4);
H = 0;
J = 1;
for i=1:count
[E(i), M(i)] = generateChain(H,J,T(i),N,iterations,MCsteps,@LCGrand);
endfor
disp(':> Zakonczono generowanie ukladu 3.');
disp([':> Czas generacji: ' num2str((t3 = cputime()) - t2)]);
figure(5);
plot(T,E)
xlabel('Temperatura');
ylabel('Energia');
title(['Energia w funkcji temperatury J = ' num2str(J) ' H = ' num2str(H) ' generator LCGrand()']);
grid;
print(['EJ' num2str(J) 'H' num2str(H) 'LCGrand.' imgExt], imgParam);
close(5);
figure(6);
plot(T,M)
xlabel('Temperatura');
ylabel('Srednia magnetyzacja');
title(['Srednia magnetyzacja w funkcji temperatury J = ' num2str(J) ' H = ' num2str(H) ' generator LCGrand()']);
grid;
print(['MJ' num2str(J) 'H' num2str(H) 'LCGrand.' imgExt], imgParam);
close(6);
H = 0;
J = -1;
for i=1:count
[E(i), M(i)] = generateChain(H,J,T(i),N,iterations,MCsteps,@LCGrand);
endfor
disp(':> Zakonczono generowanie ukladu 4.');
disp([':> Czas generacji: ' num2str((t4 = cputime()) - t3)]);
figure(7);
plot(T,E)
xlabel('Temperatura');
ylabel('Energia');
title(['Energia w funkcji temperatury J = ' num2str(J) ' H = ' num2str(H) ' generator LCGrand()']);
grid;
print(['EJ' num2str(J) 'H' num2str(H) 'LCGrand.' imgExt], imgParam);
close(7);
figure(8);
plot(T,M)
xlabel('Temperatura');
ylabel('Srednia magnetyzacja');
title(['Srednia magnetyzacja w funkcji temperatury J = ' num2str(J) ' H = ' num2str(H) ' generator LCGrand()']);
grid;
print(['MJ' num2str(J) 'H' num2str(H) 'LCGrand.' imgExt], imgParam);
close(8);
disp([':> Calkowity czas generacji: ' num2str(tt = cputime() - t0)]);
disp([':> Pojedynczy krok MC: ' num2str(tt/4/MCsteps/count)]);
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment