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

1 次查看（过去 30 天）
Jordan Lombardo
Jordan Lombardo 2018-4-12
编辑： Lee Gibson 2019-4-9
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.
  2 个评论
Lee Gibson
Lee Gibson 2019-4-9
编辑：Lee Gibson 2019-4-9
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Jim Riggs
Jim Riggs 2018-4-12
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.
Walter Roberson
Walter Roberson 2018-4-12
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.

类别

在 Help Center 和 File Exchange 中查找有关 Linear Model Identification 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by