% ### EXrandomWalk2D.m ### 11.16.14 {C. Bergevin} % Calculates mean-squared distance for a group of independent 2-D random % walkers, each moving with equal probablity a step left or right AND up or down as: % o Method 1 - simply equal probability one step up/down (U/D) and left/right (L/R) % o Method 2 - U/D and L/R steps are sampled from a uniform distribution over [-1,1] % o Method 3 - U/D and L/R steps are sampled from a Gaussian distribution % [Source (modified): Kevin Berwick (re 'Computational Physics', Giordano & Nakanishi)] clear; % ------------- N= 200; % Total # of (independent) walkers (each starts at x=0) M= 100; % Total # of steps for each walker K= 3; % # of walkers to show individual traces for [3] method= 3; % see comments above % ------------- % +++ step_number= zeros(1,M); % posAvg= zeros(1,M); % allocate array to stored (suquentially averaged) MSD step_number_array= [1:1:M]; % % +++ % NOTE: the loop is set up in such a way to average posAvg across walkers for r= 1:N x=0; y=0; % initialize positions for r'th walker positionX(r,1)= 0; positionY(r,1)= 0; % loop to go through M steps for r'th walker for nn=1:M; if method==1 if (rand<0.5), x=x+1; % conditional for left or right else x=x-1; end; if (rand<0.5), y=y+1; % conditional for up or down else y=y-1; end; elseif method==2 x= x+2*rand(1)-1; y= y+2*rand(1)-1; elseif method==3 x= x+randn(1); y= y+randn(1); end posAvg(nn)= posAvg(nn)+ (x^2+y^2); % store squared displacement (handles averging across r) positionX(r,nn+1)= x; % store displacments for each walker and step positionY(r,nn+1)= y; end; end; posAvg= posAvg/N; % Divide by number of walkers % plot MSD figure(1); plot(step_number_array,posAvg, 'k'); hold on; title('MSD for 2-D random walk'); xlabel('Step number'); ylabel('Mean-Squared Distance (x^2+y^2)'); % plot a subset of individual traces figure(2); clf; hold on; grid on; for nn=1:K shade= 1-(nn-1)/K; plot(positionX(nn,:),positionY(nn,:),'Color',[1 1 1]-shade); end xlabel('x'); ylabel('y'); title('Representative traces'); % also plot MSD (which is a circular arc in this case) rBND= sqrt(posAvg(end)); % radius MSD xBND= linspace(-rBND,rBND,100); yBND= sqrt(rBND^2-xBND.^2); plot(xBND,yBND,'r--','LineWidth',2); plot(xBND,-yBND,'r--','LineWidth',2);