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

I am unable to figure out where my error is in my code. I keep getting this message:

 Warning: Failure at t=5.469441e+00.  Unable to meet integration tolerances without reducing the step size below
the smallest value allowed (1.421085e-14) at time t.

Here is my script:

 act_thr_index = find(thr>0.25);
act_time = time(act_thr_index);
dt = length(act_thr_index)*0.0004;
act_thr = thr(act_thr_index);
g = 9.81;
mi = 43/1000;
mp = 3.12/1000;
r = 0.0132;
S = pi*r^2;
wp = mp*g;
isp = i_tot/wp;
ve = isp*g;
dv = ve*log(mi/(mi-mp));
m = linspace(mi,mi-mp,1347)';
tim = linspace(0,dt,1347)';
tspan = [0 30];
y0 = 0;
v0 = 0;
[t,y] = ode45(@odefun,tspan,[y0 v0]);

And this is my odefun:

function dydt = odefun(t,y)
%UNTITLED3 Summary of this function goes here
%   Detailed explanation goes here
dydt = zeros(2,1);
% pos = y(1);
vel = y(2);
cd = 1;
s1 = evalin('base','S');
s2 = 0.0155;
th = evalin('base','act_thr');
m = evalin('base','m');
rho = 1.225;
mi = evalin('base','mi');
mp = evalin('base','mp');
dt = evalin('base','dt');
% time = evalin('base','tim');
time = linspace(0,30,1347);
g = 9.81;
if  0 < t <= dt
  mass = interp1(time,m,t);
  thrust = interp1(time,th,t);
  drag = 0.5*rho*vel.^2.*s1*cd;
else if dt < t <= 3.7
  mass = mi-mp;
  thrust = 0;
  drag = 0.5*rho*vel.^2.*s1*cd;
  else
  mass = mi-mp;
  thrust = 0;
  drag = -0.5*rho*vel.^2.*s2*cd;
  end
dydt(1) = vel;
dydt(2) = (thrust-drag-mass*g)./mass;

I essentially used the ODEfun to compute the position and velocity of a rocket given the thrust data keeping in mind that the thrust is 0 after dt (action time) and the mass is constant after it runs out of fuel (if statement in the function).

I don't understand where to start troubleshooting.

Any help would be appreciated and i can give you more information if you need it.

[Walter Roberson]

You should not be using evalin(). See http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[Lee Gibson]

I had the opposite experience using appdesigner, doing a standard script call such as

Code

or

Code;

gave terrible and totally innacurate results, running very slowly, whereas

Code = "Code.m";

evalin('base','Code');

not only worked a lot faster, which may be part of the issue, but also gave me accurate results.

Answer 1

[Jim Riggs]

It's very hard to evaluate without all of the functions. The first thing I would check is the output of the interpolation functions to make sure that they are fairly smooth. Any large discontinuities could result in a very "stiff" system, which will be a problem for a numerical solver. This can also occur if the thrust suddenly cuts off from a large value to zero. These kinds of discontinuities can present problems to numerical solvers. (I notice that there is also a step change in the area term (s1 to s2 for T>3.7) This will also produce a discontinuity if the change is large.

If this is the case, you can break up the integral into multiple parts. Integrate up to a discontinuity (but not across) then use this final value as the initial condition for the next integral, etc.

For debugging purposes, try to artificially "smooth" out the function (i.e. set s1=s2, smooth the interpolated values, etc.) to see if this helps.

Answer 2

[Walter Roberson]

if  0 < t <= dt

means the same thing as

if ((0 < t) <= dt)

The 0 < t part returns 0 (false) or 1 (true), which you are then comparing to dt. You should be using

if 0 < t & t <= dt

However, really you should not do that at all. The ode* functions cannot handle discontinuities. You need to break the integration at each discontinuity and restart . For discontinuities based on time, the easiest way to do this is to adjust the tspan on the calls. For other discontinuities, you need to use event functions.