MATLAB: Using ode45 in a for loop using dynamic variables for particle movement and interaction

Question

Hi I am writing a small program that can create many different particles using a for loop in the x-plane that move up in the y-direction with various forces. The particles interact with each other with repelling forces and so the position of each particle with reference to each other is required in the ode45 calculation. I have looked over the code many many times but I am not sure where the fault is as the solution comes up as "NaN." The time interval is so that the movement of the particles is broken up in time so that the position of all the particles is updated in the ode45 calculation. I have used events for breaking up the time. If you need any more information as the code is not very clear, please ask! I have not shown constants in the code:

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

Any help would be really appreciated!

Thank You

Sam

EDIT

I will post a simplified version as I have taken out the segments which I know do not have any bugs:

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

Welcome to Stack Overflow, Sam! In order to make your question as answerable as possible, please reduce the code to the smallest segment that still reflects your problem. This also makes it more helpful to future readers. — nicktruesdale, Dec 05 '12 at 07:00
Hi you're right I'm sorry I will reply with a more simplified version of the code — user1878019, Dec 05 '12 at 07:10

score 1 · Accepted Answer · answered Dec 05 '12 at 07:17

1

In general, you do not want to integrate each of many equations separately, especially if the equations are coupled. ode45 and its relatives are perfectly capable of integrating a system of coupled equations.

As a simple example, you might have the following function, which defines the derivatives for three variables. The derivatives are dependent upon one another.

function xprime = my_ode(t, x)

xprime = [x(1) + 2*x(3);
          x(2) - x(1);
          3*x(3) + x(1) + x(2)];

This won't generate a very interesting system, but it shows how you can integrate three different variables in one go. Using ode45, you might type:

x0 = [0;0;1];
tspan = 0:0.1:1;
[T,X] = ode45(@my_ode, tspan, x0)

For your problem, this structure would allow you to calculate the separation between particles at each step and determine the reactions. Then, instead of integrating each and adjusting for the effects, you can make those effects a part of the equations of motion.

answered Dec 05 '12 at 07:17

nicktruesdale

815
5
12

Sorry I am slightly confused as each particle has a coupled ode of position x and y and velocity x and y. Do you mean that I should implement the position of each particle into the equations of motion? The ode45 function is in a for loop as I want to simply indicate how many particles I want, e.g. 100, and it will automatically simulate those particles and their respective motions. – user1878019 Dec 05 '12 at 07:33
If I understood right, these particles are dependent upon one another as well? If so, I would compile one long state vector, `x = [x1, y1, vx1, vy1,... xn, yn, vxn, vyn]` for `n` particles. This will allow you to fully simulate every particle and their interactions. – nicktruesdale Dec 05 '12 at 07:42
So the x in your xprime example would be 'n' long? – user1878019 Dec 05 '12 at 07:58
It would be `4*n` long to accommodate position and velocity terms. – nicktruesdale Dec 05 '12 at 17:24

MATLAB: Using ode45 in a for loop using dynamic variables for particle movement and interaction

1 Answers1