Why does ode45 output "NaN" after some time, with this time depending on the input

Hi, I am modeling a rocket launching by using ode45 to solve ODEs. I am having a strange problem where for the base problem where the rocket has 1 stage and therefore 1 input for the propellant mass, it is solved with no problem and plotted (rocket height vs. time) as shown below, but when I add another stage to the rocket and model it, ode45 outputs "NaN" starting at seemingly random timesteps.

When I change the script to account for a 2nd stage for the rocket, ode45 keeps running into instances where it can only calculate "NaN" for the solution to the odes, and this makes little sense to me because it also doesn't make sense for the derivative states to suddenly be calculated as Inf. I am able to plot the numbers it does get, and I do not understand why ode45 starts to calculate "NaN" at certain values; it almost seems random where it encounters this problem (depending on the 2 input propellant masses, of course). Shown below are the results from ode45 for the imputs of [5000;5000], [7000;7000], and [10000;0], respectively; strangely the [10000;0] case is calculated as it should be and gives a solution that I would expect. The dotted lines are the time of burnout for each stage of the rocket (so ode45 seems to have a very hard time with the second stage, but not the first). The code including the ODEs to be solved by ode45 is attached to this post as "heightcont_2stage.m". The code I was using to call the ODE function is attached to this post as "testindstudy.m".

I guess the main thing I wanted to ask was if the addition of "if" statements to govern when the equations of motion switch from the 1st stage equations to second stage equations would be what is causing this strange output, because I suspect that is what might be the problem (althought that is a guess at best). The specific code governing which stage the rocket is currently in (and what mass/acceleration to use) is shown below. There are no instances of singularities as far as I can tell either, so the reason behind the "NaN" is beyond me.

if t <= tbo1     % Burning DURING 1st stage
    mt = m0 - mdot1*t;
    vdot = (1/mt)*(Isp*g0*mdot1 - .5*rho*cd*A*x(2)^2 - (mu*mt*cos(x(5))/((RE + x(1))^2)));
elseif t > tbo1 && t <= tbo2     % Burning DURING 2nd stage
    mt = m0 - mP1 - mS1 - mdot2*t;
    vdot = (1/mt)*(Isp*g0*mdot2 - .5*rho*cd*A*x(2)^2 - (mu*mt*cos(x(5))/((RE + x(1))^2)));
elseif t > tbo2     % Coasting AFTER 2nd stage
    mt = m0 - mP1 - mS1 - mP2;
    vdot = (1/mt)*(-.5*rho*cd*A*x(2)^2 - (mu*mt*cos(x(5))/((RE + x(1))^2)));
end

I have tried solving these ODEs with several different solvers (ode45, ode15s, ode23, ode89) and they produce the same results, so it has to be something not meshing well with my code. Just wanted some guesses into what might be causing this.

Any guesses as to why this is happening are welcome, thank you.

[Torsten]

Maybe x(1) becomes -RE so that you divide by 0 ?

[Merritt]

sorry I should have added context, x(1) corresponds to the current height of the rocket off the ground and RE is the radius of the Earth (6378km). If x(1) = -RE then the rocket would be in the ground. The height values all stay above the ground before I start getting NaN, so I don't think this is the problem.

[Torsten]

clear all
close all

u_2stage = [7000;7000];
u = [10000];
tmax = 3000;
hmax = 2000;

%[heightconstraint] = heightconst(tmax,u,hmax);
[heightconstraint_2stage,xout] = heightconst_2stage(tmax,u_2stage,hmax);

Warning: Failure at t=1.383550e+02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.547474e-13) at time t.

function [constraint,xout] = heightconst_2stage(tmax,u0,hmax)
% Calculates the height of a rocket through the solution of ODEs
%
% Inputs:
%  tmax = Maximum time to run the ODE solver for (s)
%  x0 = Initial conditions for height (km), velocity (km/s), theta (rad),
%       phi (rad), and psi (rad)
%  u0 = Mass of propellant for each stage
%
% Outputs:
%  height = Height of rocket (km)

% Solve ODEs for equations of motion
tspan = [0,tmax];
odeoptions = [];
x0 = [0;0;0;0;0];
[tout,xout] = ode15s(@states,tspan,x0,odeoptions,u0);
height = xout(:,1);
maxhreached = max(height);
constraint = hmax - maxhreached;

figure(2)
plot(tout,xout(:,1),'b-')
hold on
grid on
xline(120,'k--')
xline(240,'k--')
ylim([0,5000])

for i = 1:numel(tout)
  [~,mt(i)] = states(tout(i),xout(i,:),u0);
end
figure(3)
plot(tout,mt)

% Define state space representation for use with ode45
function [output,mt] = states(t,x,u)
    % Calculates the state space representation for the objective
    % function for use in the ode45 MATLAB function
    %
    % Inputs:
    %  t = Current time step (s)
    %  x = Height (km), velocity (km/s), theta (rad), phi (rad), and psi (rad)
    %      for the previous time step
    %  u = Stage propellant masses mP1 (kg) and mP2 (kg)
    %
    % Outputs:
    %  hdot = Time derivative of height (km)
    %  vdot = Time derivative of velocity (km/s)
    %  thetadot = Time derivative of theta (rad/s)
    %  phidot = Time derivative of phi (rad/s)
    %  psidot = Time derivative of psi (rad/s)

    Isp = 353;     % Isp of propellant (s)
    epsilon = .01;     % Rocket structural ratio
    mPL = 1000;     % Payload mass (kg)     
    d = 1.65e-3;     % Radius of rocket (km)
    tbo1 = 120;     % Time to burnout for 1st stage (s)
    tbo2 = 240;     % Time to burnout for 2nd stage (s)
    cd = .3;     % Drag coefficient
    RE = 6378;     % Radius of Earth (km)
    rho0 = 1.225e9;     % Density of air at sea level (kg/km^3)
    h0 = 7.5;     % Reference height for calculation of rho (km)
    g0 = 9.81e-3;     % Gravitational acceleration at sea level (km/s^2)
    mu = 398600;     % Gravitational constant for Earth (km^3/s^2)

    mP1 = u(1,1);
    mP2 = u(2,1);
    rho = rho0*exp(-x(1)/h0);     % Air density for current height (kg/km^3)
    mS1 = mP1*(epsilon/(1 - epsilon));     % Structural mass of 1st stage (kg)
    mS2 = mP2*(epsilon/(1 - epsilon));     % Structural mass of 2nd stage (kg)
    m0 = mPL + mP1 + mS1 + mP2 + mS2;     % Initial mass of rocket (kg)
    A = (pi/4)*d^2;     % Cross-sectional area of rocket (km^2)
    mdot1 = mP1/tbo1;     % Rate of propellant use for 1st stage (kg/s)
    mdot2 = mP2/(tbo2 - tbo1);     % Rate of propellant use for 2nd stage (kg/s)
    if x(1) < 0
        x(1) = 0;
    end
    if t <= tbo1
        mt = m0 - mdot1*t;
        vdot = (1/mt)*(Isp*g0*mdot1 - .5*rho*cd*A*x(2)^2 - (mu*mt*cos(x(5))/((RE + x(1))^2)));
    elseif t > tbo1 && t <= tbo2
        mt = m0 - mP1 - mS1 - mdot2*t;
        vdot = (1/mt)*(Isp*g0*mdot2 - .5*rho*cd*A*x(2)^2 - (mu*mt*cos(x(5))/((RE + x(1))^2)));
    elseif t > tbo2
        mt = m0 - mP1 - mS1 - mP2;
        vdot = (1/mt)*(-.5*rho*cd*A*x(2)^2 - (mu*mt*cos(x(5))/((RE + x(1))^2)));
    end
    if x(2) == 0
        phidot = 0;
    else
        phidot = (mu*sin(x(5)))/(((RE + x(1))^2)*x(2));
    end
    hdot = x(2)*cos(x(5));
    thetadot = (x(2)*sin(x(5)))/(RE + x(1));
    if x(1) > .5 && x(5) < .1*(pi/180)
        psidot = .1*(pi/180);
    else
        psidot = phidot - thetadot;
    end
    output = [hdot;vdot;thetadot;phidot;psidot];
end
end

[Torsten]

It seems I found the reason for failure - you divide by mt which becomes 0 in the course of integration (see above).

[Merritt]

This seems to have been the issue! Thank you, you are a lifesaver.

Answer 1

[William Rose]

@Merritt,

There was an error when I ran testinstudy, because it calls heightconst(), which you did not provide. SO I commented out that line. It runs fine with u_2stage=[10000;0]. It runs fine with u_2stage=[5000;5000]: it doesn;t get very high, but it runs without NaNs. With u_2stage=[7000;7000], I get the result

>> testindstudy

Warning: Failure at t=2.406072e+02. Unable to meet integration tolerances without reducing

the step size below the smallest value allowed (4.547474e-13) at time t.

> In ode45 (line 353)

In heightconst_2stage (line 17)

In testindstudy (line 12)

ode45() is giving an error because it keeps trying to reduce the stepsize to get a smooth solution, and it gets to its minimum step size without finding a smooth solution. So it bails out. I have had similar problems when modeling the cardiovascular system, with one-way valves. My solution was to replace a statement like

if pressureDiff<0, 
    flow=0; 
else 
    flow=pressureDiff/resistance; 
end

with a smooth function such as

flow=pressureDiff/(resistance*(1+exp(-k*pressureDiff));

where k is large. This avoided a discontinuity in the equations, and stopped the crashing of ode45() in my simulation. I recommend you do something equivalent at the times of stage separation.

[William Rose]

@Merritt,

Elaborating on what I said above:

Replace

if t <= tbo1
    mt = m0 - mdot1*t;
    vdot = (1/mt)*(Isp*g0*mdot1 - .5*rho*cd*A*x(2)^2 - (mu*mt*cos(x(5))/((RE + x(1))^2)));
elseif t > tbo1 && t <= tbo2
    mt = m0 - mP1 - mS1 - mdot2*t;
    vdot = (1/mt)*(Isp*g0*mdot2 - .5*rho*cd*A*x(2)^2 - (mu*mt*cos(x(5))/((RE + x(1))^2)));
elseif t > tbo2
    mt = m0 - mP1 - mS1 - mP2;
    vdot = (1/mt)*(-.5*rho*cd*A*x(2)^2 - (mu*mt*cos(x(5))/((RE + x(1))^2)));
end

with

tau=0.5;  % transition duration (s)
mt = m0...
    - mdot1*t/(1+exp((t-tbo1)/tau))...      % turn off -mdot1*t when t>tbo1
    - (mP1+mS1)/(1+exp(-(t-tbo1)/tau))...   % turn on -mP1-mS1 when t>tbo1
    - mdot2*t/((1+exp(-(t-tbo1)/tau))*(1+exp((t-tbo2)/tau)))...  % turn on mdot2*t when tbo1<t<tbo2
    - mP2/(1+exp(-(t-tbo2)/tau));           % turn on -mP2 when t>tbo2
vdot = - 0.5*rho*cd*A*x(2)^2/mt - mu*cos(x(5))/(RE + x(1))^2...
       + (1/mt)*Isp*g0*mdot1/(1+exp((t-tbo1)/tau))...   % turn off +mdot1 when t>tbo1
       + (1/mt)*Isp*g0*mdot2/((1+exp(-(t-tbo1)/tau))*(1+exp((t-tbo2)/tau))); % turn on +mdot2 when tbo1<t<tbo2

Check my work, because I might have made some typos.

I don;t promise it will work, but this is the kindof think I had to do.

I am using sigmoidal funcitons that go from 0 to 1, or 1 to 0, to turn on or turn off different terms for different values of t. This gives the general idea.