Problems with ode15s

Hi. We are working on a slighty larger system in Matlab than we're used to (we are simulating interacting particles). We want to solve our differential equations for the motion of particles with matlabs ode15s. But we get the error:

??? Error using ==> feval
Undefined function or method 'Compute_ODEfunction' for input arguments of type 'double'.

Error in ==> odearguments at 110
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ==> ode15s at 228
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> particle_system>particle_system.advance_time at 79
            [~,X] = ode15s(...

Error in ==> Partikelsimulering at 64
    Particle_System.advance_time (step_time);

Our code is:

         function advance_time (Particle_System, step_time)

            Particle_System.kill_old_particles;

            time_start = Particle_System.time;
            time_end = time_start + step_time;
            tspan = [time_start, time_end];

            phase_space_state = Particle_System.get_phase_space_state';

            [~,X] = ode15s(...
                @Compute_ODEfunction, ...
                tspan, ...
                phase_space_state);

            for i = 1:6*length(Particle_System.particles)
            phase_space_state(i) = X(end,i);
            end

            Particle_System.set_phase_space_state (phase_space_state);

            Particle_System.time = time_end;

            title(Particle_System.time);

            Particle_System.advance_particles_ages (step_time);

            Particle_System.update_graphics_positions;

        end

And @Compute_ODEfunction (it's a subfunction) looks like

        function [ODEfunction] = Compute_ODEfunction ...
                (time_start, phase_space_state, Particle_System, ~, x)

            phase_space_state = phase_space_state(:)';

            Particle_System.set_phase_space_state (phase_space_state);

            Particle_System.aggregate_forces;

            NumParticles = length(Particle_System.particles);

            velocities = zeros(3*NumParticles,1);

            for i = 1:3*NumParticles
            velocities(i) = x(i+3*NumParticles);
            end

            accelerations_without_drag = Particle_System.get_particles_accelerations;
            eta = 0.0081+5.55511*time_start;
            drag_force = -Particle_System.drag*6*pi*eta*10*10^(-6);

            accelerations = zeros(3*NumParticles,1);

            for i = 1:3*NumParticles  
                accelerations(i) =  accelerations_without_drag(i) + ...
                                    drag_force*x(i+3*NumParticles)/(2*10^(-13));
            end

            ODEfunction = [velocities, accelerations]';

        end

We hope someone has an idea of what we have messed up.

Answer 1

[Walter Roberson]

To check: advance_time and Compute_ODEfunction are both defined in particle_system.m ?

I see that according to the traceback, Particle_System.advance_time (step_time) is being called, but that the next function in the traceback, particle_system>particle_system.advance_time, seems to be the advance_time function whose code you show us, which has two parameters, with the first one being Particle_System rather than step_time .

The difference in the number of parameters and the fact there is no intermediate function, tells us that Particle_System is not just a structure with a field named 'advance_time' which has been set to the handle of the function whose code we see: in order for the code to work out, Particle_System must instead be an OOP object with advance_time being one of its methods. Is that correct? If it is, then it could make a difference in the interpretation of the code, especially if Compute_ODEfunction is a private method of Particle_System . But since Compute_ODEfunction does not have Particle_System as its first parameter, Compute_ODEfunction is not a method of Particle_System at all, is it?

[Nicklas]

Yes, they are both defined in particle_system.m

And yes, that is correct - particle_system handles a lot of the different sub functions calling each other.

I didn't understood the last part? Should Compute_ODEfunction have Particle_System as its first parameter? Compute_ODEfunction is a method of Particle_System.

[Walter Roberson]

If Compute_ODEfunction is a method of a class, we are going to need to see the structure of the class definition. My recollection is that only constructors and static class methods can avoid having the object as the first argument.

For once the problem could be that, just like the message says, there is no method Compute_ODEfunction for class double -- because MATLAB only has the signature for Compute_ODEfunction of whatever class Particle_System is.

Perhaps you do indeed need to have Particle_System as the first parameter of Compute_ODEfunction, and perhaps your ode15 call should have @(x) Particle_system.Compute_ODEfunction(x)

[Nicklas]

Perfect! You were correct. It had to be called like @(x) Particle_system.Compute_ODEfunction(x). After a little more programming it worked! Thank you very much..