What kind of solver should I use for a non-continuous function?

I've been using ode45 to solve for a particular system that is easy to describe as a mass on a linear spring with no friction or other types of properties. It's a pretty simple system, but I want to implement it as being unable to become a negative length, but at the same time conserve energy (like an elastic bounce). I implemented it sort of like:

s = X(1);
v = X(2);
if s <= MinimumSpringLength
  v = abs(v)
end
~~~
Xdot = [ v, acceleration ]

My problem is, instead of actually reversing the velocity, my position hovers at MinimumSpringLength while my velocity crawls up from its initial negative state, sort of like:

X = [ 4 -4 
      2 -3
      2 -1
      2 1
      3 3 ]

These are just values I made up, but the point is that if MinimumSpringLength = 2, my position will go down to 2 and hover at 2, but has to wait for velocity to gradually accelerate to a positive. Instead, I would have preferred something like:

X = [ 4 -4 
      2 -3
      2.5 3
      4 4 ]

Is this possible to get with ode45?

Answer 1

[Are Mjaavatten]

Note that the function that you give to ode45 returns the velocity and acceleration. When you change the sign of v in this function, you correctly tell ode45 that s should increase over the next time step. The velocity change, however, is given by the acceleration, so the velocity keeps increasing, driving the position back towards minmum. You cannot effect a step change in the independent variables whithin the ode solver.

To change the velocity, you have to stop the integration when the position reaches the minimum and then restart it with the new initial velocity. You use the 'events' option to do this. Example:

 function zhang
    s0 = 5;
    v0 = 0;
    tmax = 10;
    x0 = [s0;v0];
    options = odeset('events',@reflect);
    t = 0;
    x = x0';
    while t(end) < tmax
        [t1,x1] = ode45(@spring,[t(end),tmax],x0,options);
        t = [t;t1];
        x = [x;x1];
        x0 = [x(end,1);-x(end,2)];
    end
    plot(t,x(:,1));
 end

 function dxdt = spring(t,x)
    g = -9.81;
    k = 2;
    s = x(1);
    v = x(2);
    a = g - k*(s-5);
    dxdt = [v;a];
 end

 function [value,isterminal,direction] = reflect(t,x)
    value = x(1)-1;
    isterminal = true;
    direction = -1;
 end

[Brian Zhang]

Thanks a bunch! Got it working with the event, I didn't know this existed beforehand.