## Tuesday, 9 August 2011

### Animation with MATLAB: Random walk in 2D space

This tutorial will show you how you can make a small animation to show the random walk in 2D space. Like other programs, animation is really easy. Let us see how we can do it.

Random Walk
Random walk is a walk or series of steps where each steps taken by the object is independent of previous steps. It is also known as drunken walk where the person does not know where to go and each steps taken by him is random.

We will consider the discreet case in 2D space and the steps allowed are only four as shown in the figure 1.
 Fig 1 : 2D space: Possible directions for random walk
At each time step, the object chooses a number randomly out of {1,2,3,4} and depending on that number, it will choose one of the four direction and take a step. Similar process is followed at each time step and this continues for infinite time. Random walk is the basis of many natural process like Brownian motion.

Simulation of Random Walk
Since it is an iterative process, we will need to use a loop for simulation of this. We will first simulate a single step then we will put a loop around it.

So we need to put four things.
1. assign the current position r = [x,y]
r=[0 0];
2. Get a random number B to select the direction
We will generate a number B with two binary values [b1, b2].
This two values combined will select the direction as following
$\begin{array}{|r|c|l|L|} B & & Direction & rStep \\ ----- & ----- & -------------- & ---------------- \\ [1 ,1] &-> &+X direction & [1 ,0] \\ [1 ,0 ] &-> &+Y direction & [0 ,1]\\ [0 ,1 ] &-> &-X direction &[-1 ,0]\\ [0, 0 ] &-> &-Y direction &[0 ,-1]\\ \end{array}$
B1=rand(1,2);
B=B1>0.5;
B will be 0 or 1 depending on whether B1 is less than 0.5 or not.

3. calculate new position of object with the step in the selected direction.
rnew = r + rStep;
if B==1
rnew=r+[1 0];
elseif B(1)==1
rnew=r+[0 1];
elseif B(2)==1
rnew=r+[-1 0];
else
rnew=r+[0 -1];
end
Remember if you put  B==1 as condition inside if, it will be considered true only when all the elements in B satisfies the condition.
So rnew=r+[1 0] only when B(1) and B(2) both are 1.
Otherwise if one of them is 0, it will check which one is 1 and decide on. In the end, if no conditions are met, it means that B is [0 0] and negative Y direction is chosen.

4. Update the position with this new value
%UPDATE the new position
r=rnew;
Now let us put together
r=[0 0];
B1=rand(1,2);
B=B1>0.5;

if B==1
rnew=r+[1 0];
elseif B(1)==1
rnew=r+[0 1];
elseif B(2)==1
rnew=r+[-1 0];
else
rnew=r+[0 -1];
end

%UPDATE the new position

r=rnew;
Now put a Loop around it to make it an iterative process.

r=[0 0];
for t= 0:0.1:100
B1=rand(1,2);
B=B1>0.5;

if B==1
rnew=r+[1 0];
elseif B(1)==1
rnew=r+[0 1];
elseif B(2)==1
rnew=r+[-1 0];
else
rnew=r+[0 -1];
end

%UPDATE the new position
r=rnew;
end
Here the step size is 1 and time step is 0.1 so the speed of the walk is 10 unitdistance/unittime

Animation of the walk
Only thing left is now to show the walk at each time. We will plot the step by a line which will connect the previous r to the current rnew.
plot([r(1) rnew(1)],[r(2) rnew(2)]);
Since you are adding each step on the previous trace, you need to use hold on atleast once so that previous steps are not deleted.
hold on;
plot([r(1) rnew(1)],[r(2) rnew(2)]);
Now the complete code looks like t a Loop around it to make it an iterative process.
r=[0 0];
for t= 0:0.1:100
B1=rand(1,2);
B=B1>0.5;

if B==1
rnew=r+[1 0];
elseif B(1)==1
rnew=r+[0 1];
elseif B(2)==1
rnew=r+[-1 0];
else
rnew=r+[0 -1];
end

hold on;
plot([r(1) rnew(1)],[r(2) rnew(2)]);

%UPDATE the new position
r=rnew;
end
when you run the above code, it does not give any output. The missing thing is the command 'drawnow'

Draw Now
When MATLAB sees a series of equations(commands) and a plot command in between, it calculates the graphics but it ignores the plot command at the moment to save time and when all the computation is done, then only it shows the plots. But since you want the steps to be shown at each time, we can force MATLAB to flush the plot at the moment by using drawnow command.

So the final code looks like as
r=[0 0];
for t= 0:0.1:100
B1=rand(1,2);
B=B1>0.5;

if B==1
rnew=r+[1 0];
elseif B(1)==1
rnew=r+[0 1];
elseif B(2)==1
rnew=r+[-1 0];
else
rnew=r+[0 -1];
end

hold on;
plot([r(1) rnew(1)],[r(2) rnew(2)]);

drawnow ;
%UPDATE the new position
r=rnew;
end
Run this code and you will get the following plot figure 2, You can increase time limit by changing the for loop time vector.
 Figure 2 : Random Walk
In the plot you can see that, as time passes object starts coming back to origin. It never stays at origin but it has the tendency to move around the origin and Expected translation at t= infinity is 0. This is because selection of each direction is eqaully probable since you compared with 0.5. So
$Pr[Selected Direction = +X ] = P[B(1)=1 \& B(2)=1]$ $=P(b_1 =1) \times P(b_2 =1) = P(b>0.5) \times P(b >0.5)=\frac{1}{4}$
Similarly is the other cases.
Now let us change the probabilities and see the effect.
r=[0 0];
for t= 0:0.1:100
B1=rand(1,2);
B=B1>[0.4 0.3];

if B==1
rnew=r+[1 0];
elseif B(1)==1
rnew=r+[0 1];
elseif B(2)==1
rnew=r+[-1 0];
else
rnew=r+[0 -1];
end

hold on;
plot([r(1) rnew(1)],[r(2) rnew(2)]);

drawnow ;
%UPDATE the new position
r=rnew ;
end
Now we are comparing b1 with 0.4 and b2 = 0.3. So the probability becomes
$Pr[Selected Direction = +X ] = P[B(1)=1 \& B(2)=1]$ $=P(b_1 =1) * P(b_2 =1) = P(b>0.4)*P(b >0.3)=0.6*0.7=0.42$
Similarly with other cases.

The result is shown in figure 3. You can clearly observe a tendency to move in +XY directions. This shows the motion of the object when a drifting force exists in the medium.
 Fig 3: Random walk with drift

1. Hey, thanks a lot for this project, it helped me so much! but there are an error: when you update your position, you put "rnew = r;", but this come back rnew again to zero. i put inverse way, i mean, "r = rnew;", this convert r to de current value of the coordinate. and it works for me...

2. thanks Krlos,
corrected it :)

1. Dear Gupta, Abhishek Kumar27
I need your help please. I was trying to simulate the four bar mechanism for its inverse relationship. means when the output put is considered as input. but the animation is not really same with the solution. would you please help me ? here is the code

clear;
l1=12;
l2=18;
l3=8;
l4=10;
omega1=360;
t0=(2*pi/(omega1*360));
j=1; %

for t=0:t0:((360*t0)) % the time vector
theta3(j)=(omega1*t);
A=sin(theta3(j));
B=2*(l4/l3)+cos(theta3(j));
C=(l4/l1)*cos(theta3(j))+(l4^2+l3^2+l1^2-l2^2)/(2*l3*l1);
theta1(j)=(2*atan((A+sqrt(A^2+B^2-C^2))/(B+C)));
j=j+1;
end
figure(1);
t=0:t0:(360*t0);
plot(theta3,theta1,'Erasemode','xor','linewidth',3)
title('inverse relation')
ylabel('\theta_1=f(\theta3)')
grid on;
set(gca,'Color',[1.0 1.0 1.0]);
figure(2);
plot(t,theta3)
xlabel('time')
ylabel('\theta_3')
figure(3);
plot(t,theta1)
xlabel('time')
ylabel('\theta1')

figure(2);
m=moviein(400);
i=0;
for j=1:360
i=i+1;
clf;
x(1)=0;
y(1)=0;
x(2)=l1*cos(theta1(j));
y(2)=l1*sin(theta1(j));
x(3)=-l3*cos(theta3(i));
y(3)=-l3*sin(theta3(i));
x(4)=14;
y(4)=0;
plot(x,y);
axis([-30 30 -30 30]);
title('Simulation');
grid on;
hold on;
xlabel('mm')
ylabel('mm')
plot(x(1),y(1),'o');
plot(x(2),y(2),'o');
plot(x(3),y(3),'o');
plot(x(4),y(4),'o');
m(i)=getframe;
end
movie(m);

3. That's pretty cool post! very informative and helpful, thanks for sharing this to us. We can now do this on our own.

2D Animation Philippines

4. Hi man, i need some help with writing a code for matlab. I am quite new at this and have been trying to model
a 1D random walk. I have tried to modify your code(below), but i think it is wrong. I will be extremely grateful if you could help me correct or write the code. Please

r=[0 0];
for t= 1:100
B1=rand(1,2)
if B1>0.5
rnew=r+[1 0];
else
rnew=r+[-1 0];
end
hold on;
plot(r, rnew)
drawnow ;
%UPDATE the new position
r=rnew;
end

5. hiii
i kindly need help in making a matlab animation of someone walking across a road, into a car an driving off