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

4 visualizzazioni (ultimi 30 giorni)
Merritt il 10 Nov 2023
Commentato: William Rose il 10 Nov 2023
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.
##### 5 CommentiMostra 3 commenti meno recentiNascondi 3 commenti meno recenti
Torsten il 10 Nov 2023
Modificato: Torsten il 10 Nov 2023
It seems I found the reason for failure - you divide by mt which becomes 0 in the course of integration (see above).
Merritt il 10 Nov 2023
This seems to have been the issue! Thank you, you are a lifesaver.

Accedi per commentare.

### Risposte (1)

William Rose il 10 Nov 2023
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.
##### 1 CommentoMostra -1 commenti meno recentiNascondi -1 commenti meno recenti
William Rose il 10 Nov 2023
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.

Accedi per commentare.

### Categorie

Scopri di più su Ordinary Differential Equations in Help Center e File Exchange

R2023a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by