How do I plot a graph?

36 次查看（过去 30 天）

Zaref Li 2024-8-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2148339-how-do-i-plot-a-graph

评论： Aquatris 2024-8-28

Hello everyone, I need to plot a graph according to the formulas given below. I am also sharing the image of the graph that should be. I started with the following code. I would be very happy if you could help.

reference:

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

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

Aquatris 2024-8-28

在 MATLAB Online 中打开

yes but there are infinitely many different vectors that satisfies a mean of 0 and variance of 1. So Unless you know exactly what they used to produce the graph you show, you will not be able to get the exact same result. Although we large enough simulation number, it should converge to the same solution.

Just to demonstrate, i will modify your code and run it for 100 simulations and show you the resulting first 100 x_diff and the mean:

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

% do N simulations

for n = 1:N

% 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, n) - x0) + epsilon).^5);

x_dep(i+1,n) = x_dep(i, n) + (FDEP / gamma) * dt + sqrt(2 * k_B * T / gamma) * sqrt(dt) * W(i, n); % sqrt(dt) added here

% Only diffusion (FDEP is absent)

x_diff(i+1, n) = x_diff(i, n) + sqrt(2 * k_B * T / gamma) * sqrt(dt) * W(i, n); % 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_diff * 1e6, 'b', 'LineWidth', 1.5); % y-axis scaled by 10^-6

hold on;

plot((0:num_steps-1) * dt, x_mean_diff * 1e6, 'rx', 'LineWidth', 1.5); % y-axis scaled by 10^-6

xlabel('Time (s)');

ylabel('Particle Position (μm)'); % Units updated to micrometers

title({'Particle Positions without DEP Force';'100 sim results and the mean'});

grid on;

Zaref Li 2024-8-28

编辑：Zaref Li 2024-8-28

在 MATLAB Online 中打开

First of all, thank you very much. What you said is true, I understand better now. Since W is not explicitly stated in the study, using this code will be a good choice for now. However, x_mean_dep and x_mean_diff should not give the same graph. Do you have any idea why I can't draw x_mean_dep correctly? @Aquatris

rng(0); % Reset random number generator
W = randn(num_steps, N); % Random numbers from standard normal distribution

Aquatris 2024-8-28

The Fdep value is always so small so I have no idea since I do not know the physics behind these equations.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

请先登录，再回答此问题。

类别

Sciences Physics Particle & Nuclear Physics

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by