필터 지우기
필터 지우기

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

조회 수: 17 (최근 30일)
Zaref Li
Zaref Li 2024년 8월 27일
편집: Zaref Li 2024년 10월 8일 7:31
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
Sam Chak
Sam Chak 2024년 9월 4일
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
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)

댓글을 달려면 로그인하십시오.

답변 (2개)

David Goodmanson
David Goodmanson 2024년 8월 28일
편집: David Goodmanson 2024년 8월 31일
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
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
David Goodmanson 2024년 9월 18일
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
VBBV 2024년 9월 8일
@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
Zaref Li
Zaref Li 2024년 10월 8일 7:01
@VBBV @Sam Chak Thank you very much for your comments. How can I incorporate this formula into the code?
Zaref Li
Zaref Li 2024년 10월 8일 7:16
편집: Zaref Li 2024년 10월 8일 7:31
@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.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Linear Algebra에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by