Error occurring while solving odes using ode15s

Hello,

I am trying to solve a systems of nonlinear ODE using ODE15s but I keep getting an error after about 29 seconds and I can't figure out what is causing the error.

see a screen shot of the error.

What could be the likely cause of this error.

[Walter Roberson]

Did you configure a sparse mass matrix? And does your ode routine return a single precision vector?

[Telema Harry]

@Walter Roberson thanks for the feedback.

I did not configure a sparse matrix.

I am not sure if it returns a single precision vector. How can I check that?

Summary of my experience with other solvers

ODE45

The simulation runs for the first 5000 seconds but it crashes when I run it for 10,000 seconds.

And the response of the state variables are a bit off.

ODE23

The simulation runs for 20,000 seconds and the result looks OK.

ODE113

Simulation runs for 20,000 seconds but one of the state variable produced result that is a bit off.

ODE15s

This is the preferred solver in my group and I will also like to compare the result with the other solvers.

ODE23tb

Produced similar error with ODE15s

[Walter Roberson]

Are you able to post your code for testing?

[Telema Harry]

Thanks for the inputs and your willingness to help. I wish I can easily share the code

My code is really large with several functions. I am also using weather data that is too big to share here.

[Walter Roberson]

When you invoke your ode function passing in your initial time and your initial conditions, what is class() of the result?

Also, have you configured Mass Matrix? Have you configured any special functions such as for estimation of the gradient? Is it possible that you passed in your time as single precision, or your boundary conditions as single precision?

That particular section of code, you could potentially get that error depending upon what a Mass Matrix function returned, I suspect.

[Telema Harry]

@Walter RobersonThank you once again for your feedback.

1. what is class() of the result? I am not sure. Please how can I check this?
2. have you configured Mass Matrix? I am not sure but I have included a sample of the code. The I variables are all constant while the Q_.. variables are computed at some point using the different input data.
3. Have you configured any special functions such as for estimation of the gradient? No. The functions only computes input variables required for further computation. The functions are shown in the sample code
4. Is it possible that you passed in your time as single precision? I only used time (t) once, to compute the local time. Local_time = DATE.tLST * 3600 + t;
5. your boundary conditions as single precision? I have not set boundary condition.

   options = [];
     In0 = [loc.lat; loc.lon; In.T_gas; In.T_film; In.M_gas; In.Uz; In.Z]';
      [t,output] = ode15s(@dynamics,tspan,In0, options, wind);

function dxdt = dynamics(t,y,wind)


%% VARIABLES DECLARATION
    global I DATE B
%     persistent TT
%% TRANSFER VALUES FROM y VARIABLE TO THE ACTUAL VARIABLE NAMES
      lat       = y(1);
      lon       = y(2);
      T_gas     = y(3);
      T_film    = y(4);
      M_gas     = y(5);
      Uz        = y(6);
      Z         = y(7);  
 
       [t,y']
    
    [P_air, T_air, rho_air] = atmospheremodel(Z); 
    [Uw, Vw] = interpolateWindVelocity(P_air, lon, lat, wind);
    
    Local_time = DATE.tLST * 3600 + t;
    Local_time = Local_time/3600; 
    
    [day_number, ELV] = solarangle_function(DATE, Local_time,lat,lon);
    
    dlat = M* Uw; 
    dlon = (M* Vw)/ (cos(lat*(pi/180)));
    dT_film = (Q_sun + Q_Albedo + Q_IR_planet + Q_IR_film + Q_con_ext -  Q_con_int - Q_IR_out)/(I.c_film * B.M_film);    
    alt = 25000;
    if Z > alt
        dM_gas =  -mass_flow;
    else
        dM_gas =  0;
    end
    
    
    
    dUz = (GI - B.M_gross*I.g  -  Drag_Z)/M_virtual;
    dZ = Uz;
    
      RoC =  -Uz;    
      dT_gas = Q_con_int/(I.Ye*M_gas*I.Cv) - (((I.Ye - 1)*rho_air*I.g * RoC)/(I.Ye * rho_gas*I.R_gas));

     dxdt = [dlat, dlon, dT_gas, dT_film, dM_gas, dUz, dZ]';

[Walter Roberson]

Put a breakpoint at the call to ode45. Run until there. Then command

class(tspan)
class(In0)
class(wind)
class(I)
class(B)
class(DATE)
test_output = dynamics(tspan(1), In0, wind)
class(test_output)

and tell us what results show up.

[Walter Roberson]

By the way, it is recommended that you get rid of your global variables, and that you do not use that form of ode*() call to pass extra variables. That method of passing extra variables has not been documented for over 15 years, and there are cases where it will pass the wrong information (it is incompatible with some of the newer documented options, and the documented options have priority.)

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[Telema Harry]

Thank you once again for your help. Please see below output of the code.

class(tspan) = 'double'

class(In0) = double

class (wind) = struct

class (I) = struct

class(B) = struct

class(DATE) = struct

class (test_output) = double

[Telema Harry]

@Walter Roberson Thank you so much for the recommendation on how to pass additional variables to the ode*().

I have removed all the global variables and I have re-written my function as shown below.

In0 = [loc.lat; loc.lon; In.T_gas; In.T_film; In.M_gas; In.Uz; In.Z]; 
[t, output] = Balloondynamics(tspan,In0,wind, DATE, I, B);
 
function [time, dxdt] = dynamics(tspan,y0,wind, DATE, I, B)

    [time, dxdt] = ode23(@solve_ode,tspan,y0);
    
    function dxdt = solve_ode(t,y)
      %% TRANSFER VALUES FROM y VARIABLE TO THE ACTUAL VARIABLE NAMES
      lat       = y(1);
      lon       = y(2);
      T_gas     = y(3);
      T_film    = y(4);
      M_gas     = y(5);
      Uz        = y(6);
      Z         = y(7);  
 
       [t,y']
        ....  
          
    end
    
end

I am a step closer to solving my next challenge.

So, my wind data for 1 hour is about 6.5 G and loading it takes a lot of time. Idealy, I will like to run my simulation using 1 week worth of wind data. This would mean, discarding the wind data after 1 hour simulation time (3600 seconds) and loading the next hour wind data.

e.g. load new wind data after every 3600 seconds.

I can read the wind data one after the other but I am not sure how to pass it to the ode function.

However, I was able to pass the multiple wind data to my own Runge Kutta function.

Answer 1

[Sulaymon Eshkabilov]

The computed solution might be "nan" or "inf" only. Briefly speaking the size mismatch.

[Walter Roberson]

I do not understand how this could explain the error message about trying to mix sparse and single in a * operation ?

[Telema Harry]

Thank you for your input @Sulaymon Eshkabilov

I initially thought that maybe one of the computed solution is "nan" or infinity, so I tried to debug the code just before it crashed and I discovered that all the computed solutions were all real numbers.