if statements and abs() function in variable-step ODE solvers

I was reading this post online where the person mentioned that using "if statements" and "abs()" functions can have negative repercussions in MATLAB's variable-step ODE solvers (like ODE45). According to the OP, it can significantly affect time-step (requiring too low of a time step) and give poor results when the differential equations are finally integrated. I was wondering whether this is true, and if so, why. Also, how can this problem be mitigated without resorting to fix-step solvers. I've given an example code below as to what I mean:

    function [Z,Y] = sauters(We,Re,rhos,nu_G,Uinj,Dinj,theta,ts,SMDs0,Uzs0,...
       Uts0,Vzs0,zspan,K)

    Y0 = [SMDs0;Uzs0;Uts0;Vzs0]; %Initial Conditions
    options = odeset('RelTol',1e-7,'AbsTol',1e-7); %Tolerance Levels
    [Z,Y] = ode45(@func,zspan,Y0,options);

    function DY = func(z,y)

    DY = zeros(4,1);

    %Calculate Local Droplet Reynolds Numbers
    Rez = y(1)*abs(y(2)-y(4))*Dinj*Uinj/nu_G;
    Ret = y(1)*abs(y(3))*Dinj*Uinj/nu_G;

    %Calculate Droplet Drag Coefficient
    Cdz = dragcof(Rez);
    Cdt = dragcof(Ret);

    %Calculate Total Relative Velocity and Droplet Reynolds Number
    Utot = sqrt((y(2)-y(4))^2 + y(3)^2);
    Red = y(1)*abs(Utot)*Dinj*Uinj/nu_G;

    %Calculate Derivatives
    %SMD
    if(Red > 1)
        DY(1) = -(We/8)*rhos*y(1)*(Utot*Utot/y(2))*(Cdz*(y(2)-y(4)) + ...
            Cdt*y(3)) + (We/6)*y(1)*y(1)*(y(2)*DY(2) + y(3)*DY(3)) + ...
            (We/Re)*K*(Red^0.5)*Utot*Utot/y(2);
    elseif(Red < 1)
        DY(1) = -(We/8)*rhos*y(1)*(Utot*Utot/y(2))*(Cdz*(y(2)-y(4)) + ...
        Cdt*y(3)) + (We/6)*y(1)*y(1)*(y(2)*DY(2) + y(3)*DY(3)) + ...
        (We/Re)*K*(Red)*Utot*Utot/y(2);
    end
    %Axial Droplet Velocity
    DY(2) = -(3/4)*rhos*(Cdz/y(1))*Utot*(1 - y(4)/y(2));
    %Tangential Droplet Velocity
    DY(3) = -(3/4)*rhos*(Cdt/y(1))*Utot*(y(3)/y(2));
    %Axial Gas Velocity
    DY(4) = (3/8)*((ts - ts^2)/(z^2))*(cos(theta)/(tan(theta)^2))*...
        (Cdz/y(1))*(Utot/y(4))*(1 - y(4)/y(2)) - y(4)/z;

    end

    end

Where the function "dragcof" is given by the following:

    function Cd = dragcof(Re)    
    if(Re <= 0.01)

        Cd = (0.1875) + (24.0/Re);

    elseif(Re > 0.01 && Re <= 260.0)

        Cd = (24.0/Re)*(1.0 + 0.1315*Re^(0.32 - 0.05*log10(Re)));

    else

        Cd = (24.0/Re)*(1.0 + 0.1935*Re^0.6305);
    end
    end

Answer 1

[Jan]

This has been explained in e.g. http://www.mathworks.com/matlabcentral/answers/59582#answer_72047. Matlab's ODE solvers are designed to solve functions with smooth derivatives, because the step-size control has to predict the next step based on the difference between two solutions determined with different discretizations schemas (at least for the Runge-Kutta-type solvers, but this concerns other solvers also).

The references paper of Calvo, Montijano, Randez contains a description of a method for automatic recognition of discontinuities. It is not trivial to insert it into Matlab's ODE45 and even harder for the stiff solvers. But feel free to send an enhancement request to TMW.

Problems which consider static and dynamic friction are extremely hard to simulate, because the status can switch with a very high frequency (to be exact even infinite changes are possible). Then the number of state changes determines the number of integration steps, while the dynamics of the physical body do not matter anymore. In consequence the resulting trajectories will be very sensitive and instable solutions are the standard - as you know from your personal experiences when you walked over smooth ice. Of course this is not caused by a biologically implemented Runge-Kutta intergrator, but the physical problem itself is malicious and magic is required for a solution, e.g. a neuronal network and several years of skating training.