for loop for monte carlo code

13 次查看（过去 30 天）
Elizabeth Teeter
Elizabeth Teeter 2019-12-4
评论： Elizabeth Teeter 2019-12-6
I'm coding a monte carlo simulation where a particle moves based on a vacancy next to it. I've written the step where a particle will move, but I'm unsure of how to code the for loop where in each step another random particle is selected and moves.
When i try to create the loop, i get the error that the right and left sides have a different number of elements.
Here is my code of the first step:
k=10; %total number of positions
n=k*0.5;%percentage of positions that are occupied
p=0.5; %probability of moving ys+1 if there are vacancies on both sides
x=1:k;
z=zeros(1,n,'uint32');
y=randsample(x,n); %array of occupied positions
ys=randsample(y,1); %ys is selected occupied position
find(x==ys+1);
find(x==ys-1);
Lia1=ismember(ys+1,y); %is the positive position next to ys vacant (1 is no, 0 is yes)
Lia2=ismember(ys-1,y);%is the positive position next to ys vacant (1 is no, 0 is yes)
Lia3=ismember(ys+1,x); %is ys+1 a member of x
Lia4=ismember(ys-1,x); %is ys-1 a member of x
if ys==1 && Lia2==0 %ends of array can only move inwards
yn=(y);
elseif ys==k && Lia1==0
yn=(y);
elseif Lia1==0
yn=[y(y~=ys) ys+1]; %yn is next step array where ys moves to ys+1
elseif Lia2==0
yn=[y(y~=ys) ys-1]; %yn is next step array where ys moves to ys-1
elseif Lia1==0 && Lia2==0 %if both adjacent positions are vacancies, decision which direction to move
if flip <= p
yn=[y(y~=ys) ys+1];
elseif flip >= p
yn=[y(y~=ys) ys-1];
end
elseif Lia1==1 && Lia2==1 %If neigther adjecent position is vacant, it won't move
yn=(y);
end
figure
subplot(2,1,1); plot(y,z,'o')
axis([0 10 0 10])
subplot(2,1,2); plot(yn,z,'o')
axis([0 10 0 10])
  4 个评论
显示 2更早的评论
Elizabeth Teeter
Elizabeth Teeter 2019-12-5
I'd like the loop to iterate on ys- so a selected particle that will move with each iteration.
I'd like the code to generate an array with n filled positions (x is the array and y is the positions that are filled), select one of the filled positions, move it either left or right depending on if there is a space next to it. In the next iteration it would select another filled particle (could be the one it just moved) and then move this second position.
find(x==ys+1); and find(x==ys-1); are looking at the positions next to the selected particle and then I'm using Lia to determine if that space (ys+1 or ys-1) is a member of the filled particles in the array (y).
I can generate the filled positions and move one selected filled position (ys) in my first code, but can't figure out how to code the loop to keep selecting and moving particles N times.
The z array is just so it plots on a number line, so I do want plot(y,z).
I'm not sure if this is the most efficient way to code it, but it's what I've been able to come up with.
N=2; %number of steps
k=10; %total number of positions
n=k*0.5;%percentage of positions that are occupied
p=0.5; %probability of moving ys+1 if there are vacancies on both sides
x=1:k;
z=zeros(1,n);
y=randsample(x,n); %array of occupied positions
for i=1:N
ys(i)=randsample(y(i),1); %ys is selected occupied position
find(x==ys(i)+1);
find(x==ys(i)-1);
Lia1=ismember(ys(i)+1,y); %is the positive position next to ys vacant (1 is no, 0 is yes)
Lia2=ismember(ys(i)-1,y);%is the positive position next to ys vacant (1 is no, 0 is yes)
Lia3=ismember(ys(i)+1,x); %is ys+1 a member of x
Lia4=ismember(ys(i)-1,x); %is ys-1 a member of x
if ys(i)==1 && Lia2==0 %ends of array can only move inwards
y(i)=(y);
elseif ys(i)==k && Lia1==0
y(i)=(y);
elseif Lia1==0
y(i)=[y(y~=ys(i)) ys(i)+1]; %yn is next step array where ys moves to ys+1
elseif Lia2==0
y(i)=[y(y~=ys(i)) ys(i)-1]; %yn is next step array where ys moves to ys-1
elseif Lia1==0 && Lia2==0 %if both adjacent positions are vacancies, decision which direction to move
if flip <= p
y(i)=[y(y~=ys(i)) ys(i)+1];
elseif flip >= p
y(i)=[y(y~=ys(i)) ys(i)-1];
end
elseif Lia1==1 && Lia2==1 %If neigther adjecent position is vacant, it won't move
y(i)=(y);
end
end
figure
subplot(2,1,1); plot(y,z,'o')
axis([0 10 0 10])
subplot(2,1,2); plot(y(i),z,'o')
axis([0 10 0 10])

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Ridwan Alam
Ridwan Alam 2019-12-5
This code moves one particle (randomly picked from 5) to either left or right available space. If no space available, it remains in the same space.
N=5; %number of steps
k=10; %total number of positions
n=k*0.5;%percentage of positions that are occupied
p=0.5; %probability of moving ys+1 if there are vacancies on both sides
x=1:k;
z=1:n;
ystart=randsample(x,n); %array of occupied positions
% yall holds the occupied spaces in each iteration
% ysample holds the random particle locations for each iteration
yall = ystart;
ysample = [];
y = ystart;
for i=1:N
ysample(i)=randsample(y,1); %ys is selected occupied position
right_flag = (~ismember(ysample(i)+1,y)) & ismember(ysample(i)+1,x);
left_flag = (~ismember(ysample(i)-1,y)) & ismember(ysample(i)-1,x);
if right_flag && left_flag
flip = rand;
if flip<=p
y(y==ysample(i)) = ysample(i)+1; % right-move
else
y(y==ysample(i)) = ysample(i)-1; % left-move
end
elseif right_flag
y(y==ysample(i)) = ysample(i)+1;
elseif left_flag
y(y==ysample(i)) = ysample(i)-1;
else
y = y;
end
yall = [yall;y];
end
I skipped the figure part. You will be able to handle that easily.
Hope this helps.

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Clocks and Timers 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by