Conditionals within ODE45

Question

Fawaz Zaman 2021-1-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/715918-conditionals-within-ode45

评论： Fawaz Zaman 2021-1-14

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]);
        

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Jan 2021-1-13

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 2021-1-13

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/715918-conditionals-within-ode45#answer_597308

在 MATLAB Online 中打开

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