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
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])
Elizabeth Teeter
Elizabeth Teeter 2019년 12월 6일
Yes it does, thank you!

댓글을 달려면 로그인하십시오.

채택된 답변

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 CenterFile Exchange에서 Creating and Concatenating Matrices에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by