for loop for monte carlo code

Asked by Elizabeth Teeter on 4 Dec 2019 at 21:29

Latest activity Commented on by Elizabeth Teeter on 6 Dec 2019 at 16:25

56 views (last 30 days)

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 Comments

Show 1 older comment

Ridwan Alam on 5 Dec 2019 at 16:56

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 at 17:12

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 at 16:25

Yes it does, thank you!

1 Answer

Link

Ridwan Alam 님의 답변 5 Dec 2019 17:55

채택된 답변

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.

You are now following this question

Elizabeth Teeter (view profile)