# for loop for monte carlo code

19 views (last 30 days)
Elizabeth Teeter on 4 Dec 2019
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])

Show 1 older comment
Ridwan Alam on 5 Dec 2019
what will the loop iterate on? can you show the code that threw the error?
also, just to clarify, what are the find(x==ys+1); and find(x==ys-1); doing?
and, on the plot() did you mean plot(z,y) instead of plot(y,z)?
Elizabeth Teeter on 5 Dec 2019
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])
Elizabeth Teeter on 6 Dec 2019
Yes it does, thank you!

Ridwan Alam on 5 Dec 2019
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.