How do I solve an ODE that is piecewise defined, and oscillates between the two states?

I have a mass-spring-damper model based on an example in a textbook.

It has sinusoidal input/forcing and a stiffness reduction (alpha) when the displacement is greater than 0.

My current approach

• Execute the ODE solver inside a while loop until the predetermined simulation time ( t_f ) is reached
• Use ODE event detection to trigger an event when the position results cross 0 from either side
• The "if" conditional is supposed to detect which side the spring was compressed or extended prior to reaching 0 and then switch to the other state and continue executing the loop with new initial conditions
• "odecheck" is supposed to record a 1 or a 2 depending on which function of first-order ODE's I call

Issues

• My resulting phase plot does not match up with the results from my text and I think the issue is that my code does not actually switch between the extend and compress functions
• odecheck is a string of zeros except for the very last value

Below is an image of the MSD system as well as the code i am using to accumulate the ODE solver outputs and switch between ODE functions to solve.

I've been working on this for quite a while now and I'm stumped. Any help would be greatly appreciated!

    while tout < t_f
        [T,sol,te,ye,ie] = ode15s(ode_function,[t_0 t_f],y_0,options);

        % SOURCE (L46-52): ballode
            % Accumulate output.  This could be passed as output arguments.
                nt = length(T);
                tout = [tout; T(2:nt)];
                yout = [yout; sol(2:nt,:)];
                teout = [teout; te];
                yeout = [yeout; ye];
                ieout = [ieout; ie];

        % Set new initial conditions
            y_0(1) = 0;
            y_0(2) = sol(nt,2);
        
        % Change ODE based on which side of y=0 its on
            if sol(end,1) > 0
                ode_function = @extend;
                odecheck(nt) = 1;
            elseif sol(end,1) <= 0
                ode_function = @compress;
                odecheck(nt) = 2;
            end
        % Advance tstart
            t_0 = T(nt);
    end
    

For reference: I don't think there is an issue with my ODE Functions but just in case here are the two ODEs broken down into a set of first order ODEs for the solver to interperet.

    %% ODE Functions
    function y_dot = extend(t,y)
        global m c k a w
        u = 10*sin(w*t);

        y_dot(1,1) = y(2);
        y_dot(2,1) = (u -c*y(2) - k*a*y(1))/m;
    end

    function y_dot = compress(t,y)
        global m c k w
        u = 10*sin(w*t);

        y_dot(1,1) = y(2);
        y_dot(2,1) = (u -c*y(2) - k*y(1))/m;
    end

[Jan]

Do the options contain an event function to stop the integration?

[Blake Roberts]

Currently my event function has the "isterminal" value set to 0 to continue integration, setting it to 1 seems to stop my results from concatenating within the loop.

%% y = 0 event detection
    function [position,isterminal,direction] = y0event(t,y)
        position = y(1); % The value we want to detect 0 for
        isterminal = 0;  % Continue integration (1 will halt it)
        direction = 0;   % 0 can be approached from +/- direction
    end

[Jan]

See my answer as reaktion to this comment.

Answer 1

[Jan]

The purpose of the event function is to stop the integration, such that the outer loop can restart it with the modified model. If isterminal is 0 in all cases, the event function is sleeping. Should it triggerat y(1)==0? Then:

    function [position,isterminal,direction] = y0event(t,y)
        position = y(1);
        isterminal = 1;   % !!!
        direction = 0;
    end

Answer 2

[Bjorn Gustavsson]

This modification seems to work "sensibly" OK:

Modify the ODEs into one:

function y_dot = spring_ext_comp(t,y,m,c,k,a,w,Amp)
% SPRING_EXT_COMP - 
%   Scrapped the globals - a man should have some principles in life - this
%   is my...

u = Amp*sin(w*t); 

y_dot(1,1) = y(2);
if y(1) > 0
  y_dot(2,1) = (u -c*y(2) - k*a*y(1))/m;
else
  y_dot(2,1) = (u -c*y(2) - k*y(1))/m;
end

Then I simply run for the full time-of-interest:

T = linspace(0,30,3001);
M = 1;
c = 0.1;
K = 10;
A = 1/8;
W = 10;
Amp = 0; % Just to clearly illustrate that the event-capturing properly modifies the spring-force
[T,Y,te,ye,ie] = ode15s(@(t,y) spring_ext_comp(t,y,M,c,K,A,W,Amp),T,[0.1 12],options);
sol = ode15s(@(t,y) spring_ext_comp(t,y,M,c,K,A,W,Amp),T([1 end]),[0.1 12],options);
clf
plot(T,Y(:,2))                          
hold on
plot(te,ye,'r*')                        
plot(sol.xe,sol.ye,'c.','markersize',12)
plot(sol.x,sol.y(2,:),'g')              

Which looks OK. Then you can test and check for your parameters - it also looked OK for the positions after the initial decaying natural oscillation have damped out with larger amplitudes on one side when I turned up the amplitude of the drivnig force.

HTH

[Jan]

Matlab's ODE integrators are valid, if the function to be integrated is smooth (differentiable). Otherwise the stepsize controller has an undefined behavior: Either it stop the integration with an error or the stepsize is set to such a tiny width, that the rounding errors shadow the discontinuity.

In this case, spring_ext_comp can be evaluated multiple times for one step, while some evaluations happen for y(1)<0 and some for y(1) >= 0. At least in this step the calcluated derivative is random. The effect might look small, but the scientifically correct usage would be to use an event function to stop and restart the integration at y(1)==0. Then using your single spring_ext_comp() would be fine, because the two different models are not mixed.

I've seen too many integrations in scientific publications written without asking an experienced mathematician for numerics. Integrating over discontinuities is a classical fault, which would get obvious, if the sensitivity of the result is examined: Small variations of the inputs or parameters must cause only small deviations of the final values. Otherwise the integration is instable and in consequence not a valid simulation.

[Bjorn Gustavsson]

Well, I learnt something today, something dissapointing, but still. I was assuming that if one made the effort of calculating the exact point of an event then that point would be automatically added to the solution, which would've been sufficient for this case. That was a faulty assumption.

Answer 3

[Bisesh]

%% ODE Functions

function y_dot = extend(t,y)

global m c k a w

u = 10*sin(w*t);

y_dot(1,1) = y(2);

y_dot(2,1) = (u -c*y(2) - k*a*y(1))/m;

end

function y_dot = compress(t,y)

global m c k w

u = 10*sin(w*t);

y_dot(1,1) = y(2);

y_dot(2,1) = (u -c*y(2) - k*y(1))/m;

end