3

嗨，我正在编写一个小程序，它可以使用 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
4

1 回答 1

1

通常，您不希望单独对多个方程中的每一个进行积分，尤其是在方程耦合的情况下。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)

对于您的问题，此结构将允许您计算每个步骤中粒子之间的分离并确定反应。然后，您可以将这些效果作为运动方程的一部分，而不是整合每个效果并调整效果。

于 2012-12-05T07:17:16.410 回答