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

Question

Zaref Li 2024-8-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2148124-how-can-i-scale-and-plot-the-graph-according-to-the-formula-i-use

编辑： Zaref Li 2024-10-8

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 个评论
显示 21更早的评论隐藏 21更早的评论

Sam Chak 2024-8-27

在 MATLAB Online 中打开

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 2024-8-28

在 MATLAB Online 中打开

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;

Zaref Li 2024-9-4

在 MATLAB Online 中打开

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

Sam Chak 2024-9-4

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 2024-9-8

@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 2024-8-28

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2148124-how-can-i-scale-and-plot-the-graph-according-to-the-formula-i-use#answer_1506744

编辑：David Goodmanson 2024-8-31

在 MATLAB Online 中打开

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 个评论
显示无隐藏无

Zaref Li 2024-9-4

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 2024-9-18

在 MATLAB Online 中打开

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 2024-9-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2148124-how-can-i-scale-and-plot-the-graph-according-to-the-formula-i-use#answer_1512869

在 MATLAB Online 中打开

@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 个评论
显示 3更早的评论隐藏 3更早的评论

Zaref Li 2024-10-8

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

Zaref Li 2024-10-8

编辑：Zaref Li 2024-10-8

@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.

请先登录，再进行评论。

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

23 个评论
显示 21更早的评论隐藏 21更早的评论

回答（2 个）

2 个评论
显示无隐藏无

5 个评论
显示 3更早的评论隐藏 3更早的评论

另请参阅

类别

标签

Community Treasure Hunt

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

23 个评论 显示 21更早的评论隐藏 21更早的评论

回答（2 个）

2 个评论 显示 无隐藏 无

5 个评论 显示 3更早的评论隐藏 3更早的评论

另请参阅

类别

标签

Community Treasure Hunt

23 个评论
显示 21更早的评论隐藏 21更早的评论

2 个评论
显示无隐藏无

5 个评论
显示 3更早的评论隐藏 3更早的评论