clear all;
close all;
s = 0.01; % an individual displacement of a molecule as a result of a single collision
n = input('Enter the number of molecular collisions, n = ');
t = 1 : n; % the vector of time instances
N = 500;
A = rand(N,n) - 0.5; % A is an N-by-n matrix of random numbers equally distributed in the interval [-0.5, 0.5]
A = s*sign(A);
X = zeros(N,n); % X is an N-by-n matrix of total displacements of molecules
% columns of X correspond to positions of all N molecules at a time instance
for j = 1:n
if ( j ~= 1)
X(:,j) = sum(A(:,1:j)')';
% the sum is taken for each row of A between 1 and j columns
% the result is assigned to the j column of matrix X
% X contains trajectories of all N molecules simultaneously
else
X(:,1) = A(:,1);
% no summations are performed for the first collision
end
end
Xav = sum(X)/N; % Xav is a row vector of n entries which give the mean at different time instances
XX = zeros(N,n);
for j = 1 : n
XX(:,j) = Xav(j)*ones(N,1); % the j-th column of XX has the value of Xav(j)
end
Xav2 = sum((X - XX).^2)/N; % Xav2 is a row vector of n entries which give the variance at different time instances
figure(1);
plot(Xav); % plot the mean as function of time at figure(1)
title('Evolution of the mean \mu in time');
xlabel('Time'); ylabel('\mu');
figure(2);
plot(Xav2); % plot the variance as function of time at figure(2)
title('Evolution of the variance \sigma^2 in time');
xlabel('Time'); ylabel('\sigma^2');
% computations of the coefficient of diffusion D: sigma^2 = D*n