Hitting smallest timestep allowed with several ODE solvers

Hello! I have a 6dof model in Matlab of chaff particles in a velocity field that uses an ODE solver to find the values of 12 states of the particles over time. Generally the results I get from this model are good, but occasionally the solver reaches a point in time (usually about 1.5 seconds into my 2 seconds simulation) where I get the following warning and then the solver exits:

Unable to meet integration tolerances without reducing the step size below the smallest value allowed

After scouring the internet to find people who have run into similar problems, I found the following to generally be the case and have tried each in turn without success:

1. The problem is stiff, try a different solver - I'm not certain that my problem is stiff, but I have tried ode45, ode15s, and ode23s and get the same results, along with the same warning.
2. The model is diverging/hitting a singularity, check your equations - When I plot all my states at the time of the warning, they are clearly not diverging/hitting a singularity. In fact, the warning usually occurs after all the excitement dies down and the particles I'm modelling should be simply following the velocity field at this point. The state values at the time of the warning are no farther away from the values at the previous time step than at any other point during the simulation, and certainly aren't going off to infinity.
3. The error tolerances are too tight, lower them - This potential solution I thought might fix my problem, and I went from tolerance around 1e-9 to 1e-2, but this hasn't fixed my problem. Instead, it seems like the solver just slows down to an unreasonable rate around the time it was exiting before, until it eventually gives up again after making me wait several more hours.

As far as I can tell, there seems to be no reason for the hangup. Looking at the state plots over time after it exits simply shows that the trajectories just stop, nothing blows up and there's no apparent reason for slowing down so much. My equations of motion are below:

function x_dot = eom(t,x)

Input: Time t State vector x: x(1) = u x(2) = v x(3) = w x(4) = p (roll-rate) x(5) = q (pitch-rate) x(6) = r (yaw-rate) x(7) = xi (x in inertial frame) x(8) = yi (y in inertial frame) x(9) = zi (z in inertial frame) x(10) = phi (roll) x(11) = theta (pitch) x(12) = psi (yaw)

Output: Time derivative of state vector x

Calls: atmos.m dcm.m

% Print current time
t

% Declare global variables
global g m Ixx Iyy Izz Ixy Ixz Iyz c l S use_flowfield ...
    Re_data t_data force_data lvar Fx Fy Fz alt0 vel_data

% Set chaff property values
g = 9.8;            % gravitational acceleration, m/s^2
m = 2.4494e-7;      % mass, kg (5.40e-7 lb)
Ixx = 2.5761e-16;   % roll moment of inertia, kg*m^2 (19.0e-17 slug*ft^2)
Iyy = 4.1759e-11;   % pitch moment of inertia, kg*m^2 (3.08e-11 slug*ft^2)
Izz = Iyy;          % yaw moment of inertia, kg*m^2
Ixz = 0.0;          % product moment of inertia, kg*m^2
Ixy = 0.0;          % product moment of inertia, kg*m^2
Iyz = 0.0;          % product moment of inertia, kg*m^2
c = 0.0001524;      % reference width, m (0.006 in)
l = 0.045212 + lvar;       % reference length, m (1.78 in)
S = c*l;            % reference area, m^2

Get atmospheric data

HIB = dcm(x(10),x(11),x(12));
[rho,pres,temp,a,mu] = atmos(x(9));

% Get flowfield
if use_flowfield
    ux_a = Fx(x(7),x(8),x(9)-alt0);
    uy_a = Fy(x(7),x(8),x(9)-alt0);
    uz_a = Fz(x(7),x(8),x(9)-alt0);
    if isnan(ux_a)||isnan(uy_a)||isnan(uz_a)
%         keyboard
        ux_a = 0;   uy_a = 40;  uz_a = 0;
    end
    wind = HIB'*[ux_a; uy_a; uz_a];
    u_a = wind(1);  v_a = wind(2);  w_a = wind(3);
else
    wind = HIB'*[0;0;0];
    u_a = wind(1);  v_a = wind(2);  w_a = wind(3);
end

Find forces acting on aircraft

u = x(1)-u_a;   v = x(2)-v_a;   w = x(3)-w_a;
p = x(4);       q = x(5);       r = x(6);
V = sqrt(u^2 + v^2 + w^2);
alpha = atan2(sqrt(v^2+w^2),u);
tilt = atan2(v,u);
xi = atan2(w,v);
if isnan(tilt)
    tilt = 0;
end
qbar = 0.5*rho*V^2;
qbarS = qbar*S;
Re = rho*V*c/mu;
Re_data(end+1) = Re;
t_data(end+1) = t;

[Ca,Csf,Cn,Cma,Cmn,Cmsf] = aerodata(alpha,tilt,rho,Re,V, ...
    x(4),x(5),x(6));

Cx = Ca;
Cy = Csf*sin(xi) - Cn*cos(xi);
Cz = -Cn*sin(xi) - Csf*cos(xi);
Cl = Cma;
Cm = Cmn;
Cn = Cmsf;

X = Cx*qbarS;
Y = Cy*qbarS;
Z = Cz*qbarS;
L =  Cl*qbarS*l;
M =  Cm*qbarS*c;
N =  Cn*qbarS*l;
force_data(end+1,:) = [X,Y,Z,L,M,N];

Determine the time derivatives of the state variables Based on equations 3.6-19 to 3.6-39 from Stengel

u_d = X/m - g*sin(x(11)) + x(6)*v - x(5)*w;
v_d = Y/m + g*sin(x(10))*cos(x(11)) - x(6)*u + x(4)*w;
w_d = Z/m + g*cos(x(10))*cos(x(11)) + x(5)*u - x(4)*v;

p_d = (Izz*L + Ixz*N - (Izz*(Izz - Iyy) + Ixz*Ixz)*x(5)*x(6) + ...
    Ixz*(Ixx-Iyy+Izz)*x(4)*x(5))/(Ixx*Izz - Ixz*Ixz);
q_d = (M + (Izz-Ixx)*x(4)*x(6) - Ixz*(x(4)*x(4)-x(6)*x(6)))/Iyy;
r_d = (Ixx*N + Ixz*L + (Ixx*(Ixx-Iyy)+Ixz*Ixz)*x(4)*x(5) - ...
    Ixz*(Ixx-Iyy+Izz)*x(5)*x(6))/(Ixx*Izz - Ixz*Ixz);

y = HIB * [x(1);x(2);x(3)];
x_d = y(1);
y_d = y(2);
z_d = y(3);
vel_data(end+1,:) = [x_d;y_d;z_d];

phi_d = x(4) + (x(5)*sin(x(10)) + x(6)*cos(x(10)))*tan(x(11));
theta_d = x(5)*cos(x(10)) - x(6)*sin(x(10));
psi_d = (x(5)*sin(x(10)) + x(6)*cos(x(10)))/cos(x(11));

Build state derivative vector

x_dot = [u_d; v_d; w_d; p_d; q_d; r_d; x_d; y_d; z_d; phi_d; theta_d; psi_d];

Answer 1

[Liz]

I think I figured out the problem, and its actually just a simple sign error:

When Fx, Fy, and Fz return NaN, I have a conditional to check and replace the values with (0, 40, 0) for their velocities, (ux_a, uy_a, uz_a), since it simply means the location was out of range of the Fx, Fy, and Fz functions. However, right before the particles go out of range, they have values of approximately (0, -40, 0) for their velocities. So basically, the particles were stuck at that boundary blowing back and forth instead of moving on as I intended. When I fixed the error the ODE solver was able to continue beyond that point.

[Star Strider]

Good.

Since I only have your ODE and not your ODE call or any other information, there’s no way I could have discovered that. You don’t get credit for Accepting your own Answer, so I’ll give you a vote and delete my (now irrelevant) Answer.

[John D'Errico]

Essentially, when the solver tells you it is unable to meet the integration tolerances by reducing the step size, that almost always implies your problem is stiff, in some way or another. In this case, it was stiff, due to a bug in your code.