algorithm - MATLAB:提高逻辑语句中的代码效率(Metropolis 方法)

Question

我有一个代码，旨在模拟理想气体在盒子中的运动。它基于蒙特卡罗模拟中的大都会方法。但是，我使用了一系列逻辑语句(主要是 ifs)来定义边界条件，以找到与随机选择的粒子相邻的粒子。该算法为 1x4 矩阵内的任何相邻粒子表示 1，为没有粒子的相邻点表示 0。我需要算法自动将框外的任何相邻点设置为 0，以用于框边缘上的任何粒子。有没有办法减少这些逻辑语句并仍然得到相同的结果？代码如下。

% define constants
beta = .01; %Inverse temperature
N=2000; %Duration of simulation
eps = -7;
mu=9;


for j=1:N
    a=randi(L);
    b=randi(L);
    c=randi(L);
    d=randi(L);
    % Calculate energy at positions
    if lattice(a,b)==1 && lattice(c,d)==0
        %If distribution (according to energy) suggests move,
              lattice(a,b)=0;
            lattice(c,d)=1;

         %Energy at random site/position
       E = n*(mu-eps)
            if (a~=1&& a~=L &&b~=1 && b~=L)
               adjacent= [lattice(a-1,b) lattice(a+1,b) lattice(a,b+1) lattice(a,b-1) 1];
            else if (a==1 && b==1)
                    adjacent= [0 lattice(a+1,b) lattice(a,b+1) 0 1];
                else if (a==L && b==L)
                        adjacent= [lattice(a-1,b) 0 0 lattice(a,b-1) 1];
                    else if(a==1&&b==L)
                            adjacent= [0 lattice(a+1,b) 0 lattice(a,b-1) 1];
                        else if (a==L &&b==1)
                                adjacent= [lattice(a-1,b) 0 lattice(a,b+1) 0 1];
                            else if (a==1 && b~=L && b~=1)
                                    adjacent= [0 lattice(a+1,b) lattice(a,b+1) lattice(a,b-1) 1];
                                 else if (a==L && b~=L && b~=1)
                                        adjacent= [lattice(a-1,b) 0 lattice(a,b+1) lattice(a,b-1) 1];
                                     else if (b==1&&a~=L&&a~=1)
                                             adjacent= [lattice(a-1,b) lattice(a+1,b) lattice(a,b+1) 0 1];
                                         else if (b==L&&a~=L&&a~=1)
                                              adjacent= [lattice(a-1,b) lattice(a+1,b) 0 lattice(a,b-1) 1];  
                                             end
                                         end
                                     end
                                end
                            end
                        end
                    end
                end
            end


            E1 = mu*sum(adjacent) + eps*sum(sum(adjacent.*adjacent)); 
            %This calculates the energy of the particle at its current
            %position

            if (c~=1&& c~=L &&d~=1 && d~=L)
               adjacent1= [lattice(c-1,d) lattice(c+1,d) lattice(c,d+1) lattice(c,d-1) 1];
            else if (c==1 && d==1)
                    adjacent1= [0 lattice(c+1,d) lattice(c,d+1) 0 1];
                else if (c==L && d==L)
                        adjacent1= [lattice(c-1,d) 0 0 lattice(c,d-1) 1];
                    else if(c==1&&d==L)
                            adjacent1= [0 lattice(c+1,d) 0 lattice(c,d-1) 1];
                        else if (c==L &&d==1)
                                adjacent1= [lattice(c-1,d) 0 lattice(c,d+1) 0 1];
                            else if (c==1 && d~=L && d~=1)
                                    adjacent1= [0 lattice(c+1,d) lattice(c,d+1) lattice(c,d-1) 1];
                                 else if (c==L && d~=L && d~=1)
                                        adjacent1= [lattice(c-1,d) 0 lattice(c,d+1) lattice(c,d-1) 1];
                                     else if (d==1&&c~=L&&c~=1)
                                             adjacent1= [lattice(c-1,d) lattice(c+1,d) lattice(c,d+1) 0 1];
                                         else if (d==L&&c~=L&&c~=1)
                                              adjacent1= [lattice(c-1,d) lattice(c+1,d) 0 lattice(c,d-1) 1];  
                                             end
                                         end
                                     end
                                end
                            end
                        end
                    end
                end
            end


            E2 = mu*sum(adjacent) + eps*sum(sum(adjacent1.*adjacent1)); 
            %Calculates the energy at randomly chosen position.

            dE = E2-E1; %Change in energy of the particle as it goes from the two locations.

            if rand<exp(-beta*dE)
                lattice(a,b)=0;
                lattice(c,d)=1;
            end



    end
 time(:,:,j)=lattice;


end

Answer 1

您可以取消所有这些 if else 语句。实际上，您只是在评估四个逻辑表达式并基于此选择索引。您可以使用以下内容

ab = [a-1, b; a+1, b; a, b+1; a, b-1];
abIn = [a ~= 1; a ~= L; b ~= L; b ~= 1];
adjacent = zeros(1, 5);
adjacent(abIn) = lattice(sub2ind([L, L], ab(abIn, 1), ab(abIn, 2)));
adjacent(5) = 1;

而不是第一组 if else 语句，并为第二组更改名称。这里的想法是只获取边界内的 lattice 值。四个逻辑表达式在这里用于 logical indexing在第四行。