Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
for loop ode45 with dynamic variables

Subject: for loop ode45 with dynamic variables

From: Kitech

Date: 5 Dec, 2012 02:13:11

Message: 1 of 2

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

Subject: for loop ode45 with dynamic variables

From: Nasser M. Abbasi

Date: 5 Dec, 2012 12:07:21

Message: 2 of 2

On 12/4/2012 8:13 PM, Kitech wrote:

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

Ok, so debugging is hard. But no one said programming is easy. So
you should first do as you have done, which is code inspection,
no computers. Make a nice clean print out of the code, have
a red pen in you hand, put the listing in front of you,
and read the code and make sure it is correct. Think about
the logic, see if it all flows as expected, etc...

When this does not help, the next step is to use the debugger.

Make break points, step into the code, until you see the problem show up.

It is hard work, but this is programming. debugging is part of
programming, asking someone else to debug your code for you will
not help you learn how to do it.

--Nasser

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us