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
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.