Conditionals within ODE45

Essentially, I am trying to modify the equation passed to ode45 depending on the value of y it creates. I have attempted to do this using a piecewise function, but to no avail, I get an

Error using symengine

Unable to generate code for

piecewise for use in

anonymous functions.

error.

Is there any work around this, or perhaps another method I can use to achieve this effect.

Also attached is an image that may help explain the actual physics of the problem. Thanks

function HarmonicMotion
%Default Values
    m_c = 2;        % Container Mass
    s1 = 16;        % Spring Constant 1
    s2 = 20;        % Spring Constant 2
    c = 2;          % Damping Coef
    d = 0.6;        % Gap
    t_total = 10;    % Total Integration Time
    dt = 0.01;      % Integration step
    x_0 = 2.0;      % Init Displacement
    
    tvals = [0:dt:t_total];
    x1 = nf_s1(c, m_c, s1, s2, d, dt, t_total, x_0);
    plot(tvals, x1)
    
    function [disp_vel] = nf_s1(c, m_c, s1, s2, d, dt, t_total, x_0)
        tvals = [0:dt:t_total];
        syms y(t);
        diff_x = y(t)-d; % the adj displacment for the second spring
        is_s2 = piecewise(y(t) <= d, 0, y(t) > d, 1) % turns 'on/off' the effect of the bottom spring
        
        eqn =  s1*y(t) + c*diff(y) + is_s2*s2*(diff_x) ==  -m_c * diff(y,2);
        [V] = odeToVectorField(eqn)
        M   = matlabFunction(V, 'vars', {'t', 'Y'})
        [tvals, disp_vel] = ode45(M, tvals, [x_0 0]);
        

[Jan]

Just a remark: Matlab's ODE integratore are designed to handle smooth functions only. The jump of the foirces, when the spring gets or looses contact can disturb the step size controller such that the integration stops with an error or the result is dominated by rounding errors.

Use an event function to stop the integration and restart it with the using the final value as initial value for the next interval with the changed parameters.

Answer 1

[Alan Stevens]

You could try the following simplistic approach

%Default Values
    m_c = 2;        % Container Mass
    s1 = 16;        % Spring Constant 1
    s2 = 20;        % Spring Constant 2
    c = 2;          % Damping Coef
    d = 0.6;        % Gap
    data = [m_c, s1, s2, c, d];
    t_total = 10;    % Total Integration Time
    dt = 0.01;      % Integration step
    xv0 = [2.0 0];  % Init Displacement & velocity
    
    tvals = 0:dt:t_total;
    [t, x1] = ode45(@(t,x) nf_s1(t,x,data),tvals, xv0);
    figure
    plot(t, x1(:,1),[0 t_total],[-d -d],'k--'),grid
    ylabel('displacement')
    hold on
    yyaxis right
    plot(t,x1(:,2))
    ylabel('velocity')
    xlabel('time')
    
    function [disp_vel] = nf_s1(~, xv, data)
        
       m_c = data(1); s1 = data(2); s2 = data(3); c = data(4); d = data(5);
       x = xv(1); v = xv(2);
       if x<-d
          y = x+d;
       else
           y = 0;
       end
       disp_vel = [v;
                   -(s1*x + c*v  + s2*y)/m_c]; 
      
    end

The step change in y doesn't seem large enough to adversely affect ode45 significantly here.

[Fawaz Zaman]

Hi, thanks for your solution and with neating up the code a little bit. I'm quite new to MatLab, so I do appreciate the help a lot.