我有一个代码,它是用来模拟理想气体在盒子里的运动的。它以蒙特卡罗模拟中的大都会方法为基础。然而,我使用了一系列的逻辑语句(主要是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发布于 2014-04-19 09:03:58
您可以取消所有这些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值。这四个逻辑表达式在这里用于第四行的逻辑索引。
https://stackoverflow.com/questions/23165588
复制相似问题