Matlab: if statements and abs () function in variable pitch ODE solvers

Question

I read this post on the Internet, where a person mentioned that using the "if statements" and "abs ()" functions can have negative consequences in MATLAB algorithms with a variable MATLAB step (for example, ODE45). According to the OP, it can significantly affect the time step (requiring too low a time step) and give poor results when the differential equations are finally integrated. I was wondering if this is true, and if so, why. In addition, how can this problem be mitigated without the help of solvers. I have provided the code example below, which 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 

If the dragcof function is specified as follows:

 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

This is because derivatives that are calculated using if -stations, module operations ( abs() ), or things like block step functions, delta triangles, etc., introduce gaps in the value of the solution or its derivative (s ), which leads to kinks, jumps, inflection points, etc.

This means that the ODE solution has a complete change in behavior at the appropriate times. What integrators of variable steps will do is

• define it
• recognize that they will not be able to use the information directly behind the “problem point”
• reduce the pitch and repeat from above until the problem meets the accuracy requirements

Consequently, there will be many unsuccessful steps and step size reductions near problem points, which will negatively affect the overall integration time.

However, variable step integrators will continue to get good results. Permanent step integrators are not a good way to solve this problem, because they cannot detect such problems in the first place (there is no error rating).

What you can do is simply divide the problem into several parts. If you know in advance at which points in time changes will occur, you simply start a new integration for each interval, each time using the output of the previous integration as the initial value for the next.

If you don’t know in advance where the problems will be, you can use this very nice feature in Matlab ODE solvers called event functions (see the documentation ). You allow one of the Matlab solvers to detect an event (sign change in a derivative, state change in if -statement or something else) and stop integration when such events are detected. Then start a new integration, starting from the last time and with the initial conditions of the previous integration, as before, until the end time is reached.

Thus, there will still be a small penalty in the overall execution time, since Matlab will try to pinpoint the location of the event. However, it is still much better than incorporating blind integration when it comes to runtime and accuracy of results.

Answer 2

Yes, it is true, and this is due to the fact that your decision is not smooth enough at some points.

Suppose you want to integrate. y'(t) = f(t,y) . Then what happens in f becomes integrated to become y . So if your definition of f has

• abs() , then f has a kink, and y is still smooth and 1 times differentiable
• if , then f has a jump, and y is a kink and no more differentiability.

Matlab ODE45 assumes that your solution is 5 times differentiable and tries to ensure accuracy of order 4. The nonsmooth points of your function are incorrectly interpreted as rigidity, which leads to small steps and even breakdowns.

What you can do: due to the lack of smoothness, you cannot expect high accuracy anyway. Thus, ODE23 may be the best choice. In the worst case, you must stick to first-order schemes.