Алгоритм Метрополиса в MATLAB - >> ошибка с ручками функций

Question

У меня есть фрагмент алгоритма метрополии:

mB=5.79*10^(-9); %Bohr magnetone in eV*G^-1
kB=0.86*10^(-4); %Boltzmann in eV*K^-1
%system parameters
L=60; %side square grid
L2=L*L; % total number grid position
Tstep=5; %step in temperature change (K)
Maxstep=10; %max number of steps 
nmcs=5; % cycle numberof Metropolis algorithm

magnet=NaN(1,Maxstep);%store magnetization in "monte carlo images" of sample


%Creation initial point arrangement of magnetic spins
%Outer parameters
H=100000; %Gauss
T=20; % Kelvin

%Energy alteration in spin-reverse
de =@ (i,j) (2*mB*H).*mlat(i,j);
%Metropolis probability
pmetro=@ (i,j) exp(-de(i,j)./(kB*T));

%Creation and display of initial lattice
mlat=2*round(rand(L,L))-1;

mtotal=sum(mlat(:))./L2

% Alteration of system with time
for ii=1:Maxstep
    for imc=1:nmcs
        for i=1:L
            for j=1:L
                 if pmetro(i,j)>=1
                    mlat(i,j)=-mlat(i,j);

                elseif rand<pmetro(i,j)
                    mlat(i,j)=-mlat(i,j);
                end

            end
        end
    end
    magnet(:,ii)=sum(mlat(:))./L2;
    %figure(ii);
    %pcolor(mlat);
   % shading interp;
end


m1=mean(magnet)
error=std(magnet) ./sqrt(numel(magnet))
fprintf('Temperature = %d K',T)

figure(13)
plot(magnet(1,:),'b.')
axis([0 10 0 0.5])
grid on
xlabel('i (Configuration) ')
ylabel('M/(N*mB)')

Теперь проблема в рисунке (13). Значения, которые он мне дает, равны нулю (0,05,0.02 ..). Предполагается, чтодайте мне значения около 0,3 .. В общем, график это нормально, он дает мне правильную «форму» (у него есть точки), но, как я уже сказал, около нуля.Я действительно не знаю, как поместить этот пост, чтобы его поняли. Может быть, у меня есть какая-то ошибка в матрице «магнита», я не знаю.Во всяком случае, я ни от кого не требую тщательно проверять его, я просто спрашиваю, может ли кто-нибудь помочь с быстрым взглядом.

ΕDIT --- >> Кроме того, иногда, когда я запускаю программу, она даетme:

Неопределенная функция или метод 'mlat' для входных аргументов типа 'double'.

Ошибка в ==> @ (i, j) (2 * мБ * H)). * mlat (i, j)

Ошибка в ==> @ (i, j) exp (-de (i, j) ./ (kB * T))

Ошибкав ==> мегаполисе в 39, если pmetro (i, j)> = 1

EDIT --- >>> Я нашел «ошибку». В моем коде в циклах, гдеу меня есть функция "pmetro", я заменил ее на "exp (- (2 * mB * H). * mlat (i, j) ./ (kB * T))", и программа работала просто отлично !!!Почему это не работает с вызовом «pmetro»? Как я могу это преодолеть? Есть ли проблема с дескрипторами функций в циклах?

Blockquote

Answer 1

Я настоятельно рекомендую вам попробовать написать код без использования каких-либо дескрипторов функций, пока вы действительно не познакомитесь с Matlab.

Строка de =@ (i,j) (2*mB*H).*mlat(i,j); - это то, что вызывает ваши проблемы.В Matlab, когда вы определяете дескриптор функции, который ссылается, скажем, на массив, дескриптор функции будет использовать массив таким, каким он был во время определения.Другими словами, даже если mlat изменяется внутри вашего цикла, mlat(i,j) внутри функции de всегда одинаков.На самом деле, вы даже не можете запустить этот код, если вы ранее не определили mlat в рабочей области.

Поэтому вы должны переписать основной цикл следующим образом

for iStep = 1:maxStep
   for imc = 1:mcs
      pmetro = $some function of mlat - this can be calculated using the 
                entire array as input
      %# for each element in mlat (and thus pmetro), decide whether 
      %# you have to switch the spin
      switchIdx = pmetro > 1 | pmetro < rand(size(mlat)); 
      mlat(switchIdx) = -mlat(switchIdx);
   end 
   $calculate magnetization$
end

Также обратите внимание, что тамэто команда mean, чтобы взять среднее.Не нужно суммировать, а затем делить на количество элементов.