Simulating Levy walk in MATLAB. (Not Levy Flight)

I am trying to simulate Levy Walk in 2D (Not Levy Flight). I do get the MSD and VCAF correctly, but I don't get the step length distribution correctly ( Levy Walk step lengths have heavy tails). I don't know my simulation is entirely correct. Can someone confirm my code ?

alpha=1.4;
x(1)=0;
y(1)=0;
n = 100; % 
dt =0.1; % time step
v=1; % velocity
for i=1:n   
    t= round(abs((rand()).^(-1/alpha))./dt); % time before the change in direction taken from a simple stable distribution
    theta = 2*pi*rand;  
    time(i)=t;
    for j=1:t
        x(end+1) = x(end) + v*dt*cos(theta);
        y(end+1) = y(end) + v*dt*sin(theta);
    end
end

figure(1);
plot(x, y, '-');

%% Distribution of step size
dx = diff(x);
dy = diff(y);
vxy=([dx, dy]);
bins = linspace(min(vxy), max(vxy), 75);
[count, edge]= histcounts(vxy, bins, 'Normalization', 'PDF');
pd= fitdist(transpose(vxy),'Normal');
countfit=pdf(pd, edge(1:end-1));
figure(3)
semilogy(edge(1:end-1), count, 'o', 'LineWidth',2.0); hold on;
semilogy(edge(1:end-1),countfit,'-r','LineWidth',2.0); hold off;

Answer 1

[Rahul]

Hi Alakesh,

From what I understand, you are trying to simulate a Lévy walk in 2D, but in order to achieve a better simulation result such that the step length distribution has the characteristic heavy tails of a Lévy walk, you can consider some improvements:

• Generating Levy Walk: Lévy walks are characterized by having step lengths that follow a heavy-tailed distribution, which means the step lengths L are distributed according to a power law, where the probability density function of step lengths L is proportional to L^(−(α+1)) for a Lévy exponent α, which has been done using the following line of code in the given snippet:

t=round(abs((rand()).^(-1/alpha))./dt);

• Calculating Step Lengths: Your step length distribution might be off because you are calculating the step lengths from ‘dx’ and ‘dy’ but then treating these as a single combined distribution. Since step lengths in a Levy Walk should be analyzed in terms of their magnitude, you should calculate the ‘step_lengths’ directly as:as Euclidean distances:

step_lengths = sqrt(dx.^2 + dy.^2);

Whereas your current code is doing:

vxy=([dx, dy]);  

• Plotting Histogram: You need to calculate the histogram of the step lengths and use appropriate bins to accurately capture the heavy tails, instead of calculating the histogram of ‘dx’ and ‘dy’ which are the differences between consecutive positions, since these represent the x and y components of the step, not the magnitudes.

Further, you can use ‘histcounts’ or ‘histogram’ function to plot the distribution of step-lengths and fit it to a power law, by modifying your code to plot the step length distribution:

% Plot trajectory
figure;
plot(x, y, '-');
title('Trajectory of Lévy Walk');
% Calculate step lengths
dx = diff(x);
dy = diff(y);
step_lengths = sqrt(dx.^2 + dy.^2);
% Plot histogram of step lengths
figure;
histogram(step_lengths, 'Normalization', 'pdf', 'BinEdges', logspace(log10(min(step_lengths)), log10(max(step_lengths)), 50));
title('Step Length Distribution');
xlabel('Step Length');
ylabel('Probability Density');
% Fit a power law
pd = fitdist(step_lengths', 'Lognormal');
x_values = logspace(log10(min(step_lengths)), log10(max(step_lengths)), 100);
y_fit = pdf(pd, x_values);
hold on;
semilogy(x_values, y_fit, '-r', 'LineWidth', 2.0);
hold off;

Plot the histogram of step lengths directly rather than fitting a normal distribution and use logarithmic bin edges for the histogram to better capture the heavy tails of the distribution.

For more information regarding usage of ‘logspace’ to generate logarithmically spaced vector points or log-scale plots using ‘loglog’ function in MATLAB, refer to the documentation links mentioned below:

https://www.mathworks.com/help/matlab/ref/logspace.html

https://www.mathworks.com/help/matlab/ref/loglog.html