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

6 views (last 30 days)

Show older comments

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;

##### 23 Comments

Sam Chak
on 4 Sep 2024

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
on 8 Sep 2024

### Answers (2)

David Goodmanson
on 28 Aug 2024

Edited: David Goodmanson
on 31 Aug 2024

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;

##### 2 Comments

David Goodmanson
on 18 Sep 2024

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.

VBBV
on 8 Sep 2024

@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

##### 5 Comments

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!