您好,我正在编写一个小程序,它可以在 x 平面中使用 for 循环创建许多不同的粒子,这些粒子在各种力的作用下沿 y 方向向上移动。粒子以排斥力相互作用,因此在 ode45 计算中需要每个粒子相对于彼此的位置。我已经多次查看代码,但我不确定问题出在哪里,因为解决方案显示为“NaN”。时间间隔是为了及时打散粒子的运动,以便在ode45计算中更新所有粒子的位置。我用事件来打发时间。如果您需要更多信息,因为代码不是很清楚,请询问!我没有在代码中显示常量:
time_interval = 1000; %Time in microseconds
%The number of time intervals the particle movement is divided by
number_of_intervals = 10000/time_interval;
tstart=0;
tstop= time_interval*1e-6;
num_of_particles = 3;
w0{1} = [0;0;0;0];
sols{1} = transpose(w0{1});
x_adj_stable(1) = 0;
y_adj_stable(1) = 0;
break_vect(num_of_particles) = 0;
times{1} = 0;
%Initiate particle characteristics
for particle_num = 2:num_of_particles
w0{particle_num} = [0;0;w0{particle_num - 1}(3) + 4*1e-6;0];
x_adj_stable(particle_num) = w0{particle_num}(3);
y_adj_stable(particle_num) = 0;
sols{particle_num} = transpose(w0{particle_num});
times{particle_num} = 0;
end
%Event options, Identifying the time intervals within the solution
options = odeset('Events', @events);
%The loop which connects the intervals of time
for n = 1 : number_of_intervals
for particle_num = 1:num_of_particles
%===========================================
[t_sol, y_sol] = ode45(@pm1dwoc, [tstart, tstop], w0{particle_num}, options);
%===========================================
%The number of steps within each interval
n_steps = length(t_sol);
for interval = 2: n_steps
times{particle_num} = [times{particle_num}; t_sol(interval)];
sols{particle_num} = [sols{particle_num}; y_sol(interval, :)];
end
%Setting the initial conditions of the next interval
w0{particle_num} = y_sol(n_steps,:);
y_adj(particle_num) = y_sol(n_steps, 1);
x_adj(particle_num) = y_sol(n_steps, 3);
end
tstart = tstart + time_interval*1e-6;
tstop = tstart + time_interval*1e-6;
for n = 1:num_of_particles
y_adj_stable(n) = y_adj(n);
x_adj_stable(n) = x_adj(n);
end
end
function dwdt = pm1dwoc (t,w)
y=w(1);
vy=w(2);
x=w(3);
vx=w(4);
for particle_number = 1:num_of_particles
Dist(particle_number) = sqrt((x - x_adj_stable(particle_number))^2 + (y - y_adj_stable(particle_number))^2);
if (Dist(particle_number) > 0)
x_vect(particle_number) = (x - x_adj_stable(particle_number))/Dist(particle_number);
y_vect(particle_number) = (y - y_adj_stable(particle_number))/Dist(particle_number);
end
end
dy_component = -qw*Vd/(m*G) + qw*(-qw)/(m*4*pi*epsilon0*(2*R+2*y)^2));
dx_component = 0;
for particle_number = 1:num_of_particles
if (Dist(particle_number) > 0)
dx_component = dx_component + ((qw^2)/(0.0001*m*4*pi*epsilon0*Dist(particle_number)^2))*x_vect(particle_number);
dy_component = dy_component + ((qw^2)/(0.0001*m*4*pi*epsilon0*Dist(particle_number)^2))*y_vect(particle_number);
end
end
dwdt = [vy; dy_component; vx; dx_component];
end
function [eventvalue,stopthecalc,eventdirection] = events (t,w)
t= t*1e6;
end_time = 10000;
eventvalue = t - end_time/number_of_intervals;
stopthecalc = 1;
eventdirection = 1;
for num = 2:number_of_intervals
eventvalue = [eventvalue; t - num*(end_time/number_of_intervals)];
stopthecalc = [stopthecalc; 1];
eventdirection = [eventdirection; 1];
end
end
end
非常感谢任何帮助!
谢谢
山姆
编辑
我将发布一个简化版本,因为我删除了我知道没有任何错误的部分:
num_of_particles = 3;
w0{1} = [0;0;0;0];
sols{1} = transpose(w0{1});
x_adj_stable(1) = 0;
y_adj_stable(1) = 0;
times{1} = 0;
%Initiate particle characteristics
for particle_num = 2:num_of_particles
w0{particle_num} = [0;0;w0{particle_num - 1}(3) + 4*1e-6;0];
x_adj_stable(particle_num) = w0{particle_num}(3);
y_adj_stable(particle_num) = 0;
sols{particle_num} = transpose(w0{particle_num});
times{particle_num} = 0;
end
for particle_num = 1:num_of_particles
[t_sol, y_sol] = ode45(@pm1dwoc, [tstart, tstop], w0{particle_num});
%The number of steps within each interval
n_steps = length(t_sol);
for interval = 2: n_steps
times{particle_num} = [times{particle_num}; t_sol(interval)];
sols{particle_num} = [sols{particle_num}; y_sol(interval, :)];
end
%Setting the initial conditions of the next interval
w0{particle_num} = y_sol(n_steps,:);
y_adj(particle_num) = y_sol(n_steps, 1);
x_adj(particle_num) = y_sol(n_steps, 3);
end
%//This part is written since the y_adj and x_adj variables are constantly changing so %//stable variable version is required when used in the ode45 calculations
for n = 1:num_of_particles
y_adj_stable(n) = y_adj(n);
x_adj_stable(n) = x_adj(n);
end
end
function dwdt = pm1dwoc (t,w)
y=w(1);
vy=w(2);
x=w(3);
vx=w(4);
for particle_number = 1:num_of_particles
Dist(particle_number) = sqrt((x - x_adj_stable(particle_number))^2 + (y - y_adj_stable(particle_number))^2);
if (Dist(particle_number) > 0)
x_vect(particle_number) = (x - x_adj_stable(particle_number))/Dist(particle_number);
y_vect(particle_number) = (y - y_adj_stable(particle_number))/Dist(particle_number);
end
end
dy_component = 0;
dx_component = 0;
for particle_number = 1:num_of_particles
if (Dist(particle_number) > 0)
dx_component = dx_component + (1/(Dist(particle_number))*x_vect(particle_number);
dy_component = dy_component + (1/(Dist(particle_number))*y_vect(particle_number);
end
end
dwdt = [vy; dy_component; vx; dx_component];
end
最佳答案
一般来说,您不想分别对许多方程中的每一个进行积分,尤其是如果方程是耦合的。 ode45
及其亲戚完全能够对耦合方程组进行积分。
作为一个简单的示例,您可能有以下函数,它定义了三个变量的导数。衍生品相互依赖。
function xprime = my_ode(t, x)
xprime = [x(1) + 2*x(3);
x(2) - x(1);
3*x(3) + x(1) + x(2)];
这不会生成一个非常有趣的系统,但它展示了如何一次性集成三个不同的变量。使用 ode45
,您可以键入:
x0 = [0;0;1];
tspan = 0:0.1:1;
[T,X] = ode45(@my_ode, tspan, x0)
对于您的问题,此结构允许您计算每一步中粒子之间的分离并确定 react 。然后,您可以将这些效果作为运动方程的一部分,而不是对每个效果进行积分和调整。
关于MATLAB:在 for 循环中使用 ode45,使用动态变量进行粒子运动和交互,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/13717385/