Изинг модель от Mathematica до Matlab

Question

У меня есть эта проблема от Mathematica, которая относится к модели Ising. Вот кусок кода:

J1k = Table[2 RandomInteger[] - 1, {L}, {L}];
J2k = Table[2 RandomInteger[] - 1, {L}, {L}];

energy1 := 
 Module[{ii1, ii2, jj1, jj2, kk1, kk2, dG1, ener1, ener2, ener}, 
  ener = 0.;
  Do[
   Do[
    jj1 = ii1 + 1; jj2 = ii2 + 1; kk1 = ii1 - 1; kk2 = ii2 - 1;
    If[jj1 > L, jj1 = 1]; If[jj2 > L, jj2 = 1];
    If[kk1 < 1, kk1 = L]; If[kk2 < 1, kk2 = L];
    ener1 = -J1k[[ii1, ii2]]*mlat[[ii1, ii2]]*mlat[[jj1, ii2]] - 
      J1k[[kk1, ii2]]*mlat[[ii1, ii2]]*mlat[[kk1, ii2]];
    ener2 = -J2k[[ii1, ii2]]*mlat[[ii1, ii2]]*mlat[[ii1, jj2]] - 
      J2k[[ii1, kk2]]*mlat[[ii1, ii2]]*mlat[[ii1, kk2]];
    ener = ener1 + ener2 - 2*mlat[[ii1, ii2]] *\[Mu]Bk,
    {ii1, 1, L}],
   {ii2, 1, L}];
  ener = ener/2.]

energy = energy1

Вот что я сделал:

J1k=2*randint(L,L)-1;
J2k=2*randint(L,L)-1;
I have created a function:
function Energy1 =energy1()
%spin interaction with the neighbors
ener=0;
L=16;
J1k=2*randint(L,L)-1;
J2k=2*randint(L,L)-1;
for ii2=1:L
for ii1=1:L
    jj1=ii1+1;jj2=ii2+1;kk1=ii1-1;kk2=ii2-1;
    if (jj1>L & jj1==1)
    end
        if (jj2>L & jj2==1)
        end
            if (kk1<1 & kk1==L)
            end
                if (kk2<1 & kk2==L)
                end
                    ener1=-J1k(ii1,ii2).*mlat(ii1,ii2).*mlat(jj1,ii2)-J1k(kk1,ii2).*mlat(ii1,ii2).*mlat(kk1,ii2);
                    ener2=-J2k(ii1,ii2).*mlat(ii1,ii2).*mlat(ii1,jj2)-J2k(ii1,kk2).*mlat(ii1,ii2).*mlat(ii1,kk2);
                    ener=ener1+ener2-2.*mlat(ii1,ii2).*mBk;


end                
end
ener=ener ./2;

end

energy =@ energy1

Проблема в том, что когда я называю energy "energy = @ energy1", значение не возвращается. Предполагается, что оно возвращает значение "1" или "-1".

Кроме того, онопродолжает (в Mathematica):

T = 5.2; dT = 0.1; Umean = {};
energy = energy1;
Do[energyseries = {}; T = T - dT;
  Do[i1 = RandomInteger[{1, L}]; i2 = RandomInteger[{1, L}];
   j1 = i1 + 1; j2 = i2 + 1; k1 = i1 - 1; k2 = i2 - 1;
   If[j1 > L, j1 = 1]; If[j2 > L, j2 = 1];
   If[k1 < 1, k1 = L]; If[k2 < 1, k2 = L];
   dG1 = -J1k[[i1, i2]]*mlat[[j1, i2]] - 
     J1k[[k1, i2]]*mlat[[k1, i2]];
   dG2 = -J2k[[i1, i2]]*mlat[[i1, j2]] - 
     J2k[[i1, k2]]*mlat[[i1, k2]];
   dE = 2.*mlat[[i1, i2]] (dG1 + dG2 + \[Mu]Bk);
   W = N[Exp[-dE/T]];
   If[W < 1 && W > Random[] || dE <= 0, 
    mlat[[i1, i2]] = -mlat[[i1, i2]]; index = 1, index = 0];
   energy = energy + dE*index // N;
   AppendTo[energyseries, energy],
   {i, 1, Maxstep}];
  PrependTo[Umean, {T, Mean[energyseries]/L2}],
  {jT, 1, 51}];

И я сделал:

T=5.2;
dT=0.1;
Umean=zeros();
energy=@energy1
for jT=1:51
    energyseries=zeros();
    T=T-dT;
    for i=1:Maxstep
    i1=randint(1,L);
    i2=randint(1,L);
    j1=i1+1;
    j2=i2+1;
    k1=i1-1;
    k2=i2-1;
    if (j1>L & j1==1)
    end
        if (j2>L & j2==1)
        end
            if (k1<1 & k1==L)
            end
                if (k2<1 & k2==L)
                end
                    dG1=-J1k(i1,i2).*mlat(j1,i2)-J1k(k1,i2).*mlat(k1,i2);
                    dG2=-J2k(i1,i2).*mlat(i1,j2)-J2k(i1,k2).*mlat(i1,k2);
                    dE=2.*mlat(i1,i2).*(dG1+dG2+mBk);
                    W=exp(-dE./T);
                    if (W<1 & W>rand() | dE<=0)
                        mlat(i1,i2)=-mlat(i1,i2);
                        index=1;
                        index=0;
                    end
                    energy=energy+dE.*index;
                    energyseries(:)=[energy]
    end
    Umean(:)=[T mean(energyseries(:))./L2]
end

, что дает мне сообщение "Индексы индекса должны быть либо действительными положительными целыми числами или логическими.Ошибка в ==> исходит при 58 dG1 = -J1k (i1, i2). * Mlat (j1, i2) -J1k (k1, i2). * Mlat (k1, i2); "Есть идеи?

Answer 1

Я реализовал модель Ising в Matlab несколько лет назад. Возможно, мой код будет вам полезен; это доступно в этой заметке: Исследование Монте-Карло модели Изинга (PDF) . Код очень короткий и Matlabesque и начинается на странице 6.

В качестве примера идиоматического программирования на Matlab предположим, что у вас есть матрица grid, которая кодирует вашу сетку спинов как +1 (вверх) или -1 (вниз). Затем вы можете вычислить среднюю намагниченность соседей каждого сайта без использования циклов for, выполнив:

% Calculate the number of neighbors of each cell
neighbors = circshift(grid, [ 0 1]) + ...
            circshift(grid, [ 0 -1]) + ...
            circshift(grid, [ 1 0]) + ...
            circshift(grid, [-1 0]);

Затем вы можете вычислить изменение энергии за счет переворота любых спинов:

% Calculate the change in energy of flipping a spin
DeltaE = 2 * J * (grid .* neighbors);

Это утверждение работает на каждом узле сетки одновременно; такая «векторизация» (вычисление сразу над целыми массивами) является ключом к эффективному коду Matlab.

Answer 2

Для вашей первой проблемы вы не установили значение для выхода "Energy1".Просто сделайте

Energy1 = ener;

в конце функции.

Я думаю, что вашей второй проблемой является использование операторов if.«Если» в Matlab не работает, как «если» в Mathematica.В Matlab это работает так:

if (expression==true)   
     execute statement
end

Итак, что вы хотите сделать, это

if (j1>L)
     j1 = 1;
end

и т. Д.

РЕДАКТИРОВАТЬ, чтобы суммировать обсуждение в комментариях: Используйте«rand (L)» для определения i1 и i2 вместо «randint»