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;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;Β 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;Β graph, graphics, mathematics, matlab, loopΒ MATLAB Answers β New Questions
β
