# Generating a discrete RW on a unit circle and calculate its long-time distribution

6 views (last 30 days)
Evangelos Mitsokapas on 7 Oct 2017
Hi, I am new to MATLAB and I would like to generate a random walk of a particle on a ring (circle of radius r). For each iteration, the particle should move to the "right" on the circle by 1, with a probability p=1/2 and to the "left" on the circle by -1 with a probability q=1-p=1/2. I can generate a random walk in MATLAB like in the following code, but I am not sure how to confine it on a circle. After this, I would like to calculate the long-time distribution of the particle numerically.
%Random Walk in 2D
close all;
clear all;
nsteps=10;
x(1)=0;
y(1)=0;
for i=1:nsteps
theta=2*pi*rand()
r=1;
dx=r*cos(theta);
dy=r*sin(theta);
x(i+1)=x(i)+dx;
y(i+1)=(x(i)/tan(theta))+dy;
end
plot(x,y,'r-');

John D'Errico on 7 Oct 2017
Ok. New to MATLAB or not, I'll give you credit for having made a credible effort.
So what are you doing wrong here? You are looking at it from the wrong direction, trying to solve the problem in a bass-ackwards way. What does that mean?
This is NOT a random walk in the (x,y) plane! This is a random walk in polar angle theta, and the (x,y) coordinates are computed parameters as a function of that polar angle theta.
What you are trying to do is make it a walk in the plane, then constrain the points to lie on a circle. That is coming at the problem the wrong way around.
So just generate the angle theta as a sequence, then compute x and y. This ensures that the radius is ALWAYS 1, so the points ALWAYS lie exactly on the unit circle.
r = 1;
ntheta = 100000;
theta0 = 0;
dtheta = rand(1,ntheta-1)*2 - 1;
theta = cumsum([theta0,dtheta]);
x = r*cos(theta);
y = r*sin(theta);
plot(x,y,'.')
grid on
axis equal
axis square
I'll let you do the plots and the computations.
As far as computing the long term distribution, consider the difference between these two calls to hist:
hist(theta,100)
hist(mod(theta,2*pi),100)
Which is appropriate for your question?
Finally, could you do this using a loop? Of course. But you need to learn to use MATLAB, as it was designed to be used. To solve the problem using a loop, just create theta incrementally. Then at the end, compute x and y.
Image Analyst on 8 Oct 2017
The vertical axis is the frequency (count, number of occurrences). The x axis is the value. So it says how many times each number occurs. You can use the newer function histogram() instead of hist() if you have a recent release. It has several normalization options, like 'count' (default), 'probability', 'countdensity', 'pdf', 'cumcount', and 'cdf'.
Evangelos Mitsokapas on 8 Oct 2017
@Image Analyst: Thank you for the answer! Now it becomes more clear what has been going on behind the histogram. Would you be so kind as to provide an example of normalization of this histogram using the histogram() command? Many thanks.

Image Analyst on 7 Oct 2017
Edited: Image Analyst on 7 Oct 2017
Your first step is right on the edge of the circle! So you need to have a different step length than the size of the outer, limiting circle. Then you need to computer y(i+1) differently (just add dy to it) and have a while loop to try again if the point would land outside the outer circle.
numberOfSteps = 2000;
maxIterations = 1000;
x(1) = 0;
y(1) = 0;
stepLength = 1;
for k = 1 : numberOfSteps
again = 1;
loopCount = 1;
while again && loopCount < maxIterations
theta=2*pi*rand();
dx = stepLength * cos(theta);
dy = stepLength * sin(theta);
x(k+1) = x(k) + dx;
y(k+1) = y(k) + dy;
if x(k+1)^2 + y(k+1)^2 <= r^2
again = 0;
end
loopCount = loopCount + 1;
end
end
plot(x,y,'r.-', 'MarkerSize', 18);
grid on;
xlabel('X');
ylabel('Y');
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
ax.FontSize = 30;
ax.FontWeight = 'Bold';
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
##### 3 CommentsShow 1 older commentHide 1 older comment
Image Analyst on 7 Oct 2017
Oh, sorry - I misunderstood.
Evangelos Mitsokapas on 7 Oct 2017
It's totally fine! Again, many thanks for your time :)