How do I plot a graph?

Hello everyone, I need to plot a graph according to the formulas given below. I am also sharing the image of the graph that should be. I started with the following code. I would be very happy if you could help.

reference:

% Define parameters
epsilon0 = 8.85e-12;    % F/m (Permittivity of free space)
epsilon_m = 79;         % F/m (Relative permittivity of the medium)
CM = 0.5;               % Clausius-Mossotti factor
k_B = 1.38e-23;         % J/K (Boltzmann constant)
R = 1e-6;               % m (Particle radius)
gamma = 1.88e-8;        % kg/s (Friction coefficient)
q = 1e-14;              % C (Charge of the particle)
dt = 1e-3;              % s (Time step)
T = 300;                % K (Room temperature)
x0 = 10e-6;              % Initial position (slightly adjusted from zero)
N = 100000;             % Number of simulations
num_steps = 1000;       % Number of steps (simulation for 1 second)
epsilon = 1e-9;         % Small offset to avoid division by zero
k = 1 / (4 * pi * epsilon_m); % Constant for force coefficient

% Generate random numbers
rng(0); % Reset random number generator
W = randn(num_steps, N); % Random numbers from standard normal distribution

% Define position vectors (with and without DEP force)
x_dep = zeros(num_steps, N);
x_diff = zeros(num_steps, N);
x_dep(1, :) = x0;
x_diff(1, :) = x0;

% Perform iterations using the Euler-Maruyama method
for i = 1:num_steps-1
    % With DEP force (FDEP is present)
    FDEP = 4 * pi * R^3 * epsilon0 * epsilon_m * CM * (-2 * k^2 * q^2) ./ ((abs(x_dep(i, :) - x0) + epsilon).^5);
    x_dep(i+1, :) = x_dep(i, :) + (FDEP / gamma) * dt + sqrt(2 * k_B * T / gamma) * sqrt(dt) * W(i, :);  % sqrt(dt) added here
    
    % Only diffusion (FDEP is absent)
    x_diff(i+1, :) = x_diff(i, :) + sqrt(2 * k_B * T / gamma) * sqrt(dt) * W(i, :);  % sqrt(dt) added here
end

% Calculate means
x_mean_dep = mean(x_dep, 2);
x_mean_diff = mean(x_diff, 2);

% Plot results (y-axis scaled by 10^-6)
figure;
plot((0:num_steps-1) * dt, x_mean_dep * 1e6, 'b', 'LineWidth', 1.5); % y-axis scaled by 10^-6
hold on;
plot((0:num_steps-1) * dt, x_mean_diff * 1e6, 'r--', 'LineWidth', 1.5); % y-axis scaled by 10^-6
xlabel('Time (s)');
ylabel('Particle Position (μm)'); % Units updated to micrometers
title('Particle Positions with and without DEP Force');
legend('DEP Force Present', 'Only Diffusion');
grid on;

[Malay Agarwal]

Hi @Zaref Li,

What is the initial value of and ?

[Zaref Li]

I updated the code. I cannot draw the Fdep graphic. @Malay Agarwal

[Aquatris]

The W vector you use has significant impact on the results so in order to get the exact same graph, you need to supply exactly the same W. Does the paper mention how they created the W vector?

[Zaref Li]

expressed as "Wn​ is a random variable with a Gaussian distribution having a mean of 0 and a variance of 1." @Aquatris

[Aquatris]

yes but there are infinitely many different vectors that satisfies a mean of 0 and variance of 1. So Unless you know exactly what they used to produce the graph you show, you will not be able to get the exact same result. Although we large enough simulation number, it should converge to the same solution.

Just to demonstrate, i will modify your code and run it for 100 simulations and show you the resulting first 100 x_diff and the mean:

% Define parameters
epsilon0 = 8.85e-12;    % F/m (Permittivity of free space)
epsilon_m = 79;         % F/m (Relative permittivity of the medium)
CM = 0.5;               % Clausius-Mossotti factor
k_B = 1.38e-23;         % J/K (Boltzmann constant)
R = 1e-6;               % m (Particle radius)
gamma = 1.88e-8;        % kg/s (Friction coefficient)
q = 1e-14;              % C (Charge of the particle)
dt = 1e-3;              % s (Time step)
T = 300;                % K (Room temperature)
x0 = 10e-6;              % Initial position (slightly adjusted from zero)
N = 100;             % Number of simulations
num_steps = 1000;       % Number of steps (simulation for 1 second)
epsilon = 1e-9;         % Small offset to avoid division by zero
k = 1 / (4 * pi * epsilon_m); % Constant for force coefficient

% Generate random numbers
rng(0); % Reset random number generator


W = randn(num_steps,N) ; % Random numbers from standard normal distribution
% Define position vectors (with and without DEP force)
x_dep = zeros(num_steps, N);
x_diff = zeros(num_steps, N);
x_dep(1, :) = x0;
x_diff(1, :) = x0;

% do N simulations
for n = 1:N
% Perform iterations using the Euler-Maruyama method
for i = 1:num_steps-1
    % With DEP force (FDEP is present)
    FDEP = 4 * pi * R^3 * epsilon0 * epsilon_m * CM * (-2 * k^2 * q^2) ./ ((abs(x_dep(i, n) - x0) + epsilon).^5);
    x_dep(i+1,n) = x_dep(i, n) + (FDEP / gamma) * dt + sqrt(2 * k_B * T / gamma) * sqrt(dt) * W(i, n);  % sqrt(dt) added here
    
    % Only diffusion (FDEP is absent)
    x_diff(i+1, n) = x_diff(i, n) + sqrt(2 * k_B * T / gamma) * sqrt(dt) * W(i, n);  % sqrt(dt) added here
end
end
% Calculate means
x_mean_dep = mean(x_dep, 2);
x_mean_diff = mean(x_diff, 2);
% Plot results (y-axis scaled by 10^-6)
figure;
plot((0:num_steps-1) * dt, x_diff * 1e6, 'b', 'LineWidth', 1.5); % y-axis scaled by 10^-6
hold on;
plot((0:num_steps-1) * dt, x_mean_diff  * 1e6, 'rx', 'LineWidth', 1.5); % y-axis scaled by 10^-6
xlabel('Time (s)');
ylabel('Particle Position (μm)'); % Units updated to micrometers
title({'Particle Positions without DEP Force';'100 sim results and the mean'});
grid on;

[Zaref Li]

First of all, thank you very much. What you said is true, I understand better now. Since W is not explicitly stated in the study, using this code will be a good choice for now. However, x_mean_dep and x_mean_diff should not give the same graph. Do you have any idea why I can't draw x_mean_dep correctly? @Aquatris

rng(0); % Reset random number generator
W = randn(num_steps, N); % Random numbers from standard normal distribution

[Aquatris]

The Fdep value is always so small so I have no idea since I do not know the physics behind these equations.