19 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])

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.

Opportunities for recent engineering grads.

Apply Today
## 4 Comments

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/494899-for-loop-for-monte-carlo-code#comment_774579

⋮## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/494899-for-loop-for-monte-carlo-code#comment_774579

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/494899-for-loop-for-monte-carlo-code#comment_774864

⋮## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/494899-for-loop-for-monte-carlo-code#comment_774864

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/494899-for-loop-for-monte-carlo-code#comment_774871

⋮## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/494899-for-loop-for-monte-carlo-code#comment_774871

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/494899-for-loop-for-monte-carlo-code#comment_775299

⋮## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/494899-for-loop-for-monte-carlo-code#comment_775299

Sign in to comment.