Sorry- I just edited the question!
for loop for monte carlo code
13 views (last 30 days)
Show older comments
Elizabeth Teeter
on 4 Dec 2019
Commented: Elizabeth Teeter
on 6 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])
4 Comments
Accepted Answer
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.
0 Comments
More Answers (0)
See Also
Categories
Find more on Creating and Concatenating Matrices in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!