USER
Modify the MatLab code below to answer the following prompt.
Adjust pw by gradually increasing it from 0 to near 1. As you do this, repeat number 7 in the previous problem for different values of pw. Determine the appropriate way to translate number of steps to time by modifying the relationship: t = n𝜏 (the new formula will be a function of pw ). Now, consider two kinds of media. Media A fills the pipe for x ≥ 0, while media B fills the pipe for x < 0. Suppose the diffusion constant for media A is 10^-7 cm^2 / sec and the diffusion constant for media B is 10^-8 cm^2 / sec. By using different pw for media A and B, answer the following questions.
1) Visualize the random walk trajectories
2) In the x-direction, calculate the mean squared displacement as a function of number of steps.
3) In the y-direction, calculate the mean squared displacement as a function of number of steps.
4) For 1,000,000 walks, plot the histogram of locations of the random walker in the x-direction after 100, 1000, 10000, 100000, 1000000 steps.
5) For 1,000,000 walks, plot the histogram of locations of the random walker in the y-direction after 100, 1000, 10000, 100000, 1000000 steps.
%% RANDOM WALKTS, BROWNIAN MOTION, AND DIFFUSION
% ————————— parameters
width = 5; % width of pipe in mm
a = 0.01; % lattice constant in mm
pw = 0; % probability of waiting
num_steps = 1000000; % # steps for histograms
num_walks = 1000000; % # random walks for histograms
% ————————— calculate lattice size based on width
lattice_size = round(width/a);
% ————————— initialize random walker position
x = 0;
y = lattice_size/2;
% ————————— initialize mean squared displacement variables
msd_x = zeros(1,num_steps);
msd_y = zeros(1,num_steps);
% ————————— loop over the # of steps
for step = 2:num_steps+1
% choice: does walker wait or make a step?
if rand<=pw
% walker waits
continue;
else
% walker makes a step
% calculate the step direction
dx = randi([-1,1]);
dy = randi([-1,1]);
% update walker position
x = x+dx;
y = y+dy;
% reflect at the y-direction boundary
if y<1
y = 2-y;
elseif y>lattice_size
y = 2lattice_size - y - 1;
end
% calculate mean squared displacement in each direction
msd_x(step) = (xa)^2;
msd_y(step) = ((y-(lattice_size/2))a)^2;
end
end
% —————————— 1) visualize the random walk trajectories
figure;
plot((0:num_steps)a, (msd_ya), ‘b.’)
title(‘Random Walk Trajectories’);
xlabel(‘x (mm)’);
ylabel(‘y (mm)’);
xlim([0,num_stepsa]);
% ————————— 2) mean squared displacement in x-direction
figure;
n_steps = 1:num_steps;
msd_x = cumsum(msd_x) ./ (1:num_steps+1);
plot(n_steps(1:num_steps-1)a, msd_x(2:num_steps)a^2, ‘b’)
title(‘Mean Squared Displacement in the x-Direction’);
xlabel(‘Number of Steps’);
ylabel(‘Mean Squared Displacement’);
xlim([0,num_stepsa]);
% ————————— 3) mean squared displacement in y-direction
figure;
msd_y = cumsum(msd_y) ./ (1:num_steps+1);
plot(n_steps(1:num_steps)a, msd_y(1:num_steps)a^2, ‘b’)
title(‘Mean Squared Displacement in the y-Direction’);
xlabel(‘Number of Steps’);
ylabel(‘Mean Squared Displacement’);
xlim([0, num_stepsa]);
% ————————— 4) histogram of locations in x-direction
figure;
for num_steps = [100,1000,10000,100000,1000000]
x_hist = zeros(1,lattice_size+1);
for walk = 1:num_walks
x=0;
for step = 1:num_steps
dx = randi([-1,1]);
x = x+dx;
if x<0
x=0;
elseif x>lattice_size
x = lattice_size;
end
end
x_hist(x+1) = x_hist(x+1) + 1;
end
subplot(2,3,find([100,1000,10000,100000,1000000] == num_steps));
bar(0:lattice_size, x_hist / num_walks, ‘b’);
title(['Histogram of x-locations after ‘, num2str(num_steps), ’ steps’]);
xlabel(‘x’);
ylabel(‘Probability’);
end
% ————————— 5) histogram of locations in y-direction
figure;
for num_steps = [100,1000,10000,100000,1000000]
y_hist = zeros(1,lattice_size+1);
for walk = 1:num_walks
x = 0;
y = lattice_size/2;
for step = 1:num_steps
dx = randi([-1,1]);
dy = randi([-1,1]);
x = x+dx;
y = y+dy;
if y<1
y = 2-y;
elseif y>lattice_size
y = 2lattice_size - y - 1;
end
end
y_hist(y+1) = y_hist(y+1) + 1;
end
subplot(2,3,find([100,1000,10000,100000,1000000] == num_steps));
bar(0:lattice_size, y_hist / num_walks, ‘b’);
title(['Histogram of y-locations after ‘, num2str(num_steps), ’ steps’]);
xlabel(‘y’);
ylabel(‘Probability’);
end
% ————————— 6) calculate avg time per step for given diffusion constant
D_desired = 1e-7; % desired diffusion constant in cm^2/sec
t = D_desired / (a^2); % avg time per step in sec
% substitude the values, t = (10^-7) / (.01^2) = 10^-3
% so the average time per step (t) should be 10^-3 seconds
% ————————— 7) calculate diffusion constant in the x-direction
figure;
diffusion_x = gradient(msd_xa^2) / (2t);
plot(n_stepsa, diffusion_x);
title(‘Diffusion Constant in the x-Direction’);
xlabel(‘Time (mm)’);
ylabel(‘Diffusion Constant’);
xlim([0, num_stepsa]);
ASSISTANT
%% RANDOM WALKTS, BROWNIAN MOTION, AND DIFFUSION
% ————————— parameters
width = 5; % width of pipe in mm
a = 0.01; % lattice constant in mm
pw_A = 0; % probability of waiting for media A
pw_B = 0; % probability of waiting for media B
num_steps = 1000000; % # steps for histograms
num_walks = 1000000; % # random walks for histograms
D_A = 10^-7; % diffusion constant for media A (cm^2/sec)
D_B = 10^-8; % diffusion constant for media B (cm^2/sec)
% ————————— calculate lattice size based on width
lattice_size = round(width/a);
% ————————— initialize random walker position
x = 0;
y = lattice_size/2;
% ————————— initialize mean squared displacement variables
msd_x = zeros(1,num_steps);
msd_y = zeros(1,num_steps);
% ————————— loop over the # of steps
for step = 2:num_steps+1
% choice: does walker wait or make a step?
if rand<=pw_A
% walker waits
continue;
else
% walker makes a step
% calculate the step direction
dx = randi([-1,1]);
dy = randi([-1,1]);
% update walker position
x = x+dx;
y = y+dy;
% reflect at the y-direction boundary for media A
if y<1 && x >= 0
y = 2-y;
elseif y>lattice_size && x >= 0
y = 2*lattice_size - y - 1;
end
% reflect at the y-direction boundary for media B
if y<1 && x < 0
y = 2-y;
elseif y>lattice_size && x < 0
y = 2*lattice_size - y - 1;
end
% calculate mean squared displacement in each direction
if x >= 0
msd_x(step) = (x*a)^2;
end
msd_y(step) = ((y-(lattice_size/2))*a)^2;
end
end
% —————————— 1) visualize the random walk trajectories
figure;
plot((0:num_steps)*a, (msd_y*a), 'b.')
title('Random Walk Trajectories');
xlabel('x (mm)');
ylabel('y (mm)');
xlim([0,num_steps*a]);
% ————————— 2) mean squared displacement in x-direction
figure;
n_steps = 1:num_steps;
msd_x_cum = cumsum(msd_x) ./ (1:num_steps+1);
plot(n_steps(1:num_steps-1)*a, msd_x_cum(2:num_steps)*a^2, 'b')
title('Mean Squared Displacement in the x-Direction');
xlabel('Number of Steps');
ylabel('Mean Squared Displacement');
xlim([0,num_steps*a]);
% ————————— 3) mean squared displacement in y-direction
figure;
msd_y_cum = cumsum(msd_y) ./ (1:num_steps+1);
plot(n_steps(1:num_steps)*a, msd_y_cum(1:num_steps)*a^2, 'b')
title('Mean Squared Displacement in the y-Direction');
xlabel('Number of Steps');
ylabel('Mean Squared Displacement');
xlim([0, num_steps*a]);
% ————————— 4) histogram of locations in x-direction for media A
figure;
for num_steps = [100,1000,10000,100000,1000000]
x_hist_A = zeros(1,lattice_size+1);
for walk = 1:num_walks
x=0;
for step = 1:num_steps
dx = randi([-1,1]);
x = x+dx;
if x<0
x=0;
elseif x>lattice_size
x = lattice_size;
end
end
x_hist_A(x+1) = x_hist_A(x+1) + 1;
end
subplot(2,3,find([100,1000,10000,100000,1000000] == num_steps));
bar(0:lattice_size, x_hist_A / num_walks, 'b');
title(['Histogram of x-locations in media A after ', num2str(num_steps), ' steps']);
xlabel('x');
ylabel('Probability');
end
% ————————— 4) histogram of locations in x-direction for media B
figure;
for num_steps = [100,1000,10000,100000,1000000]
x_hist_B = zeros(1,lattice_size+1);
for walk = 1:num_walks
x=0;
for step = 1:num_steps
dx = randi([-1,1]);
x = x+dx;
if x>0 || x < 0
continue;
end
end
x_hist_B(x+1) = x_hist_B(x+1) + 1;
end
subplot(2,3,find([100,1000,10000,100000,1000000] == num_steps));
bar(0:lattice_size, x_hist_B / num_walks, 'b');
title(['Histogram of x-locations in media B after ', num2str(num_steps), ' steps']);
xlabel('x');
ylabel('Probability');
end
% ————————— 5) histogram of locations in y-direction for media A
figure;
for num_steps = [100,1000,10000,100000,1000000]
y_hist_A = zeros(1,lattice_size+1);
for walk = 1:num_walks
x = 0;
y = lattice_size/2;
for step = 1:num_steps
dx = randi([-1,1]);
dy = randi([-1,1]);
x = x+dx;
y = y+dy;
if y<1
y = 2-y;
elseif y>lattice_size
y = 2*lattice_size - y - 1;
end
end
y_hist_A(y+1) = y_hist_A(y+1) + 1;
end
subplot(2,3,find([100,1000,10000,100000,1000000] == num_steps));
bar(0:lattice_size, y_hist_A / num_walks, 'b');
title(['Histogram of y-locations in media A after ', num2str(num_steps), ' steps']);
xlabel('y');
ylabel('Probability');
end
% ————————— 5) histogram of locations in y-direction for media B
figure;
for num_steps = [100,1000,10000,100000,1000000]
y_hist_B = zeros(1,lattice_size+1);
for walk = 1:num_walks
x = 0;
y = lattice_size/2;
for step = 1:num_steps
dx = randi([-1,1]);
dy = randi([-1,1]);
x = x+dx;
y = y+dy;
if y<1 && x < 0
y = 2-y;
elseif y>lattice_size && x < 0
y = 2*lattice_size - y - 1;
end
end
y_hist_B(y+1) = y_hist_B(y+1) + 1;
end
subplot(2,3,find([100,1000,10000,100000,1000000] == num_steps));
bar(0:lattice_size, y_hist_B / num_walks, 'b');
title(['Histogram of y-locations in media B after ', num2str(num_steps), ' steps']);
xlabel('y');
ylabel('Probability');
end
% ————————— 6) calculate avg time per step for media A
t_A = D_A / (a^2); % avg time per step for media A in sec
% ————————— 7) calculate diffusion constant in the x-direction for media A
figure;
diffusion_x_A = gradient(msd_x_cum*a^2) / (2*t_A);
plot(n_steps*a, diffusion_x_A);
title('Diffusion Constant in the x-Direction for media A');
xlabel('Time (mm)');
ylabel('Diffusion Constant');
xlim([0, num_steps*a]);
% ————————— 6) calculate avg time per step for media B
t_B = D_B / (a^2); % avg time per step for media B in sec
% ————————— 7) calculate diffusion constant in the x-direction for media B
figure;
diffusion_x_B = gradient(msd_x_cum*a^2) / (2*t_B);
plot(n_steps*a, diffusion_x_B);
title('Diffusion Constant in the x-Direction for media B');
xlabel('Time (mm)');
ylabel('Diffusion Constant');
xlim([0, num_steps*a]);