Math Forum Discussions

Kitech

Posts: 2
Registered: 12/4/12

for loop ode45 with dynamic variables
Posted: Dec 4, 2012 9:13 PM

Plain Text

Hi I am writing a small program that can create many different particles that move up in the y-direction spread out in the x-direction. 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." 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