How can I scale and plot the graph according to the formula I use?

Hello everyone,

I am trying to plot only the diffusion and FDEP diffusion using the formulas below. However, the FDEP diffusion plot is incorrect. I am also attaching the reference graph that it should match.

formulas:

reference graph:

my graph:

% 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;

[Torsten]

The initial position in the reference graph is not 0 ...

[Zaref Li]

How can I modify my graph to match the reference graph? @Torsten

[Torsten]

How can I modify my graph to match the reference graph?

I don't know. I just saw that the curves in the reference graph start at 1e-5 whereas yours start at 0. This should be easy to modify in your code.

Further, k is undefined in your code.

[Zaref Li]

I defined k in the code. It may be easy but this is the part where I'm stuck :) I'd be glad if you could help. @Torsten

[Voss]

Your code is missing a sqrt(dt) in the calculation of x_dep(i+1, :) (and I guess the calculation of x_diff(i+1, :) as well).

But the real problem is that when x_dep(i, :) is equal to x0 (as it is initially) FDEP becomes -Inf because of dividing by zero (i.e., x_dep(i, :) - x0 is zero).

[Mathieu NOE]

@Zaref Li

not only does your graph start at the correct value, but also your y range is way different from the reference picture

error in the code or constants ?

[Zaref Li]

I tried to make the changes you suggested in the code and updated it, but unfortunately, there was no change. If you can make any adjustments, please feel free to let me know.@Voss

As you can see, the graph follows the same pattern as the reference graph after x exceeds 0.4. I don't understand how this is happening. I am just trying to apply the formula. @Mathieu NOE

[Voss]

@Zaref Li: I tried a few things, but I didn't find a way around the divide-by-0 problem, at least not one that gave a plot like the reference plot.

[Sam Chak]

Hi @Zaref Li

I would suggest double-checking your initial value x0 again. Additionally, it appears that either the formula for the parameter k is incorrect, or the formula presented in the manuscript has been incorrectly typed. If the manuscript contains an error in such a simple formula, then there is a possibility that Equations (18), (19), and (20) may also contain errors.

When the noise term W is disabled, I would be interested in understanding the expected response computed by the Euler–Maruyama numerical method.

% 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 = 1e-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
W = zeros(num_steps, N);

% 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
figure;
plot((0:num_steps-1) * dt, x_mean_dep, 'b', 'LineWidth', 1.5);
hold on;
plot((0:num_steps-1) * dt, x_mean_diff, 'r--', 'LineWidth', 1.5);
xlabel('Time (s)');
ylabel('Particle Position (μm)'); % Units updated to micrometers
title('Particle Positions');
legend('DEP Force Present', 'Only Diffusion', 'location', 'east');
grid on;

[Zaref Li]

I couldn't see where you made changes in the code you provided. Could you please explain? In your code, it starts from 10 as we wanted. @Sam Chak

I made such a change in the graph. It resembles the reference graph a bit. However, it plots the graph under the effect of the Fdep force incorrectly. Can you see why? @Torsten @Voss @Mathieu NOE @Sam Chak

% Define parameters
epsilon0 = 8.85e-12;    % F/m (Permittivity of free space)
epsilon_m = 79;         % F/m (Relative permittivity of 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 set to 10 micrometers (10x10^-6 m)
N = 100000;             % Number of simulations
num_steps = 1000;       % Number of steps (simulation for 1 second)
epsilon = 1e-9;         % Small offset to prevent 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) included
    
    % Only diffusion (FDEP is absent)
    x_diff(i+1, :) = x_diff(i, :) + sqrt(2 * k_B * T / gamma) * sqrt(dt) * W(i, :);  % sqrt(dt) included
end

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

% Determine the index corresponding to 0.4 seconds
start_index = round(0.4 / dt); % 0.4 seconds to index conversion

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

[Sam Chak]

@Zaref Li, I only disabled the randomness:

W = zeros(num_steps, N);

[Zaref Li]

Hello again, everyone. I learned that starting the initial position from 0 is not an issue. So, as far as I understand, there is currently no problem with plotting the diffusion graph. However, it is plotting the Fdep graph incorrectly. Can you help me find the problem here? I am updating the question. @Torsten @Voss @Mathieu NOE @Sam Chak

[Ronit]

Hi @Zaref Li, could you please provide the link to the paper you are referring to?

[Zaref Li]

Hi @Ronit Unfortunately, this is not taken from a paper. These are our own study notes. I am sharing them in detail.

[Sam Chak]

Hi @Zaref Li, Please show the original continuous-time differential equation of the eletric charge system. Probably Equations (11) to (16).

[Zaref Li]

Will these be useful to you? @Sam Chak

[Sam Chak]

@Zaref Li, Because Eqs. (15) and (16) are decoupled, are you able solve these two simple equations first?

[Zaref Li]

I didn't quite understand what you meant. Do I need to apply these formulas in the code? (16) is already solved with (17) @Sam Chak

[Sam Chak]

Hi @Zaref Li

Put other equations aside. From the image you supplied, are you able to use ode45() solve Eqs. (15) and (16)?

[Zaref Li]

I got help from ChatGPT for this. @Sam Chak

For the 1st Equation:

function dxdt = first_ode(t, x, m, gamma, F_ext, zeta)
dxdt = zeros(2,1);
dxdt(1) = x(2); % dx/dt = v
dxdt(2) = (-gamma/m) * x(2) + (1/m) * F_ext(x(1)) + (1/m) * zeta(t);
end

% Parameters
m = 1; % Mass (example value)
gamma = 0.1; % Friction coefficient (example value)
F_ext = @(x) -x; % Example external force function (e.g., F_ext(x) = -k*x)
zeta = @(t) 0; % Random force (example: zero)

% Initial conditions [x(0), v(0)]
x0 = [0; 0]; % x(0) = 0, v(0) = 0

% Time span
tspan = [0 10];

% Solution
[t, x] = ode45(@(t, x) first_ode(t, x, m, gamma, F_ext, zeta), tspan, x0);

% Plot results
plot(t, x(:,1))
xlabel('Time (s)')
ylabel('Position x(t)')
title('Solution of First Equation')
grid on

For the 2nd Equation:

γdx(t)dt=Fext(x)+ζ(t)\gamma \frac{dx(t)}{dt} = F_{\text{ext}}(x) + \zeta(t)γdtdx(t)​=Fext​(x)+ζ(t)

This is already a first-order equation.

function dxdt = second_ode(t, x, gamma, F_ext, zeta)
dxdt = (1/gamma) * (F_ext(x) + zeta(t));
end

% Parameters
gamma = 0.1; % Friction coefficient (example value)
F_ext = @(x) -x; % Example external force function (e.g., F_ext(x) = -k*x)
zeta = @(t) 0; % Random force (example: zero)

% Initial condition
x0 = 0; % x(0) = 0

% Time span
tspan = [0 10];

% Solution
[t, x] = ode45(@(t, x) second_ode(t, x, gamma, F_ext, zeta), tspan, x0);

% Plot results
plot(t, x)
xlabel('Time (s)')
ylabel('Position x(t)')
title('Solution of Second Equation')
grid on

[Torsten]

ode45 or any other deterministic integrator cannot be used to integrate stochastic differential equations.

[Sam Chak]

Hi @Zaref Li

Thank you for your update. I had thought you were satisfied with @David Goodmanson's solution following several days of silence.

However, I am not familiar with how you are solving the Brownian-like motion using the algorithm in Eq. (20). I was attempting to point out that when Eq. (16) is substituted into Eq. (15) and simplified, I am also unsure about how to proceed, as your ChatGPT-generated code suggests that you may be unfamiliar with solving ODEs. I believe this is perfectly acceptable for beginners. I also began to take MATLAB more seriously, exploring various tools during the COVID period.

Perhaps you could try using a third-party Euler–Maruyama algorithm directly from the File Exchange before working on your code.

[VBBV]

@Zaref Li, In your reference graph, the legend representation for blue and red curves seems incorrect. The blue curve i presume is ONLY with diffusion (which is shown as with FDEP) since it follows a straight line pattern, and the red curve is sloping downwards which is actually the FDEP (shown as diffusion only)

Answer 1

[David Goodmanson]

Hi Zaref,

One thing for sure about eq.(19) in the attachment is that eps0 can't be in the numerator. eps_m is dimensionless (although the comment in the your code says the units are F/m). Dimensionally, with eps0 in the denominator then FDEP looks like

x^3 q^2 / (eps0 x^5) = q^2 / (eps0 x^2) 

which has correct units of force. Given that the initial expression is incorrect it is not totally clear where the factors of eps_m and 4pi should go in the entire expression including k. 4pi usually accompanies eps0, so for FDEF I will use

FDEP = C*(R^3*/(4*pi*eps0))*CM*(-2*q^2) / abs(x- x1)^5                (1)

where C contains dimensionless factors of 4pi and eps_m, whatever they may be. In the code below I used C = 1. Wih a suitable value for other constants (see text below) the result is

< I temporarily dropped the 1e6 in the plot and forgot to put it back. The y axis units are meters so the plot y axis description is incorrect>.

When you average the displacement over 1e5 runs, the random motion compared to that of a single particle is going to be reduced by a factor of 1/sqrt(1e5). The drift velocity due to FDEP is not random at all, and there is really no need to average it over runs since it never changes. So for an average over a lot of runs, the importance of diffusion vs. drift will be underemphasized compared to doing just a single run.

Looking at the reference graph, the drift velocity is about -4e-9 and if you want to reproduce that, then using x=0 in (1) and

v = FDEP / gamma             

you will arrive at x1 = 1e-4, approximately. Of course given v and with x=0 in (1) you can always backsolve for x1 exactly. And if your C is different from 1 you can do the same thing to get a different x1.

The FDEP force is attractive and the drift velocity in the reference plot is negative, so the charged particle goes to the left of the origin, at -1e-4. In 1 sec the drift to the left is negligilble compared to that distance.

The code below uses a starting distance of 0 and drops some other startup constants since there are no issues at the origin. And I temporarily dropped the 1e6 in the plot and forgot to put it back, so the y axis units are meters and the plot description is incorrect.

% 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 = 0;                 % Initial position (slightly adjusted from zero)
x1 = -1e-4;      % position of charged particle
N = 100000;             % Number of simulations
num_steps = 1000;       % Number of steps (simulation for 1 second)
k = 1 / (4 * pi * epsilon_m); % Constant for force coefficient
C = 1;
% suitable drift velocity
Cdriftv = C*(R^3/(4*pi*epsilon0))*CM*(-2*q^2)/x1^5/gamma
% 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 = C*(R^3/(4*pi*epsilon0))*CM*(-2*q^2)./(abs(x_dep(i, :) - x1).^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(1);
plot((0:num_steps-1) * dt, x_mean_dep * 1e0, 'b', 'LineWidth', 1.5); % y-axis scaled by 10^-6
hold on;
plot((0:num_steps-1) * dt, x_mean_diff * 1e0, 'r--', 'LineWidth', 1.5); % y-axis scaled by 10^-6
hold off
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;

[Zaref Li]

First of all, thank you very much for adjusting the code, but I couldn't fully understand what you meant. The reference graph can be created using only the formulas I sent, and I want to create the same one. As far as I understand, the drift velocity is not needed. Also, the graph you plotted works inversely compared to the reference graph. Could you please reevaluate? @David Goodmanson

[David Goodmanson]

Hi Zaref,

I didn't see your comment until today. Could you mention the spots that are not clear? The equations I used are the same as yours. It's just a question of what to use for the constants.

As a couple of people have pointed out, it appears that the graph that is posted in the original question has the traces swapped. (I had not noticed that before). In that plot, the blue curve is relatively flat and the red curve has a noticeable negative slope. That means that the blue curve is diffusion only, and the red curve has a contribution from FDEP.

THe FDEP contribution is (I used C=1)

x(j+1)-x(j) = C*(R^3*/(4*pi*eps0))*CM*(-2*q^2) / abs(x- x1)^5 dt      

This term is a negative factor times dt, meaning tha there is a negative drift velocity.

The particle starts at the origin and in the times of interest it barely moves at all compared to x1. So the expression can be taken as

x(j+1)-x(j) = C*(R^3*/(4*pi*eps0))*CM*(-2*q^2) / abs(x1)^5 dt

and in the times of interest the drift velocity is (approximately) constant.

Answer 2

[VBBV]

@Zaref Li you need to take the differential for noise term while solving the equation i.e. dW instead of W

% 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+1, :)-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+1, :)-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

[Zaref Li]

but it doesn't look like this reference chart? @VBBV

[VBBV]

@Zaref Li, According to the algorithm given in snapshot of your question, info from this link https://en.wikipedia.org/wiki/Euler%E2%80%93Maruyama_method shows that for noise term, differential term needs to be used. From your code, this has not been considered/.

[Sam Chak]

Nice catch! I revisited the Wiki page and compared the formula with the one in the snapshot.

[Zaref Li]

@VBBV @Sam Chak Thank you very much for your comments. How can I incorporate this formula into the code?

[Zaref Li]

@VBBV @Sam Chak I indicated it as (W(i+1, :) - W(i, :)) in the code, isn't this correct? But we still can't match the reference graph. It's showing the opposite behavior.