Problem with returning to Main function from subfunction

Hi everybody,

I'm running into a problem when I want to execute a code that is schematic like the following (See below). It is supposed to produce the orbital decay of space objects. I use a ODE integrator like Runge-Kutte Fehleberg to get the rates of change. When a certain value is reached (like altitude is 0) I want to abort the integration and return the values to the main function for plotting etc. The integration and orbit decay seems to work ok but when the abort criteria is reached I get the following error:

-------------------------------------------------------------------

_Error in SimpleDecay>rates (line 135) dfdt = zeros(7,1);

Output argument "dydt" (and maybe others) not assigned during call to "C:\Users\...\SimpleDecay.m>SimpleDecay/rates".

Error in rkf45 (line 79) f(:,i) = feval(ode_function, t_inner, y_inner);

Error in SimpleDecay (line 86) [t,f] = rkf45(@rates, tspan, f0);

_--------------------------------------------------------------

I'm not sure whether I placed the abort 'if' correctly or what the problem really is so every help is highly appreciated!

Thanks in advance, David

The code is schematic like this:

function main

%Some definitions etc.

%Define timespan of integration & initial value vector

% Calling the integrator

[t,f] = rkf45(@rates, tspan, f0);

%return the solutions

a = f(:,1);

%etc.

plotting %as a subfunction call

return

function  dydt = rates(ti,y)

dfdt = zeros(7,1);

initial vector definition

a = y(1);

%etc.

% Define some needed values

% Return when value is reached

if a < Re

      return

end

%Calculate rates

dadt = ...

% etc.

% Load derivatives into output vector

dydt(1) = dadt;

end % end rates

end % end main

[Andrew Newell]

This ODE solver is from MATLAB Central. Have you tried contacting the author? Also, have you considered using the MATLAB solvers like ode45?

[David]

Thanks for the fast reply. Not so far, and for my understanding, but I could be wrong, the ODE is not the problem because is works as expected when I keep the simulation time short enough to come not near the abort condition. What I don't see is the reason why the program does not "simple" abort, return the values to the main function and plots it when I simulate until the criterion is reached.

Answer 1

[Walter Roberson]

You must assign values to any output variable used by the calling function if you are going to use "return" (explicitly or implicitly). This is basic MATLAB, not dependent on the fact that you are using an ode solver.

Early stopping for the MathWorks-provided ode* solvers is handled by passing an options structure (usually created with odeset) that includes an Events option that provides a function handle to test for a stop.

Examining the code for the rkf45 contribution, I see that it has no provision for early stopping except for this: returning a very large value from the user function will cause early termination. The easiest way to get the necessary large value would be to return infinity.

[David]

Thank you that was the right hint!

Answer 2

[David]

Thanks a lot for that hint with odeset Walter. I figured how to do that and it works under certain circumstances. General an event location property should help with that. Here is the simplified code I wrote with the help of some helpful internet resources:

http://www.mathworks.com/help/techdoc/ref/ode23.html

http://www.mathworks.com/moler/odes.pdf

http://www.mathworks.com/help/techdoc/ref/odeset.html#f92-1016858

The code:

function test
t0 = 0;
tf =  5e7;
tspan = [t0,tf];
a0 = 6578.135;
f0 = a0;
e = 0;
mu = 398600.79964;
Re = 6378.135;

A = 38.5;
m = 135e3;
cD = 2.214;
BC = cD*(A/m);

opts=odeset('Events',@events);

[t,f] = ode45(@rates, tspan, f0, opts);
fprintf('ODE finish\n');

a = f(:,1);

plot(t, a-Re,'-x')
xlabel('Time (?)')
ylabel('Semi-major axis in [a-Re] (km)')
axis([-inf, inf, -inf, inf])
grid
return

function  dydt = rates(ti,y)
    dfdt = zeros(1,1);
    a = y(1);
        n = sqrt(mu/a^3);
        h = a - Re;
H = 8.4; %km
rho0 = 1.225; %kg/m^3
rho_pe = rho0*exp(-h/H);
        dadt = -BC*rho_pe*1000*(a^2)*(n/(2*pi))*1*(1+(3/4)*(e^2)*(1));
        fprintf('Rates calculated\n');
        fprintf('a\th\n');
        fprintf('%10g %10g \r\n',a,a-Re);

      dydt(1) = dadt;
      fprintf('End rates\n');
  end

      function [value,isterminal,direction] = events(t,y)
  % Locate the time when height passes through zero in a decreasing
  % direction
  % and stop integration.
  value = y(1)-(Re+190); % detect a-Re = 0
  fprintf('Abort?\n');
  isterminal = 1; % stop the integration
  direction = []; % negative direction
      end
  end

It works well with the value of "190" in the "value" line (6th from the end). But when I e.g. change it to 100 the plot stops below 170. And the command window output shows that the last values are far below 100 but shouldn't the ODE been already stopped? This causes also some issues with additional calls in the sub functions "rates" for e.g. the atmospheric density when I use an other external function for that. My understanding was that ODE stops right after the value 100 is reached (or a value near that). Also a marker ("Abort?") seems not been called after every operation of the ODE. The question is now how can stop the ODE right after reaching the abort value or what is my mistake that the code doesn't do it? Hope it is clear what I mean and somebody has an answer or hint! Thanks in advance!

Answer 3

[David]

Finally I found a solution. It is actually a numerical issue. I think the integrator oversees the abort value and is anyway to fast in the lower atmosphere. This is solved by increasing the 'RelTol' default value by e.g. 1e-12 with a odeset option.

Thank you Walter and Andrew for there contribution! I really appreciate that.

The final code of my example is (with some other changes):

function test
t0 = 0;
tf =  5e7;
tspan = [t0,tf];
a0 = 6578.135;
f0 = a0;
e = 0;
mu = 398600.79964;
Re = 6378.135;

A = 38.5;
m = 135e3;
cD = 2.214;
BC = cD*(A/m);
refine = 1; %How many intermediate values
opts=odeset('Events',@events,'RelTol',1e-12,'Refine',refine); 

%Defaults: RelTol = .001 and AbsTol = .000001

[t,f] = ode45(@rates, tspan, f0, opts);
  if ~ishold
      hold on
  end
   nt = length(t);
  opts = odeset(opts,'InitialStep',t(nt)-t(nt-refine),...
          'MaxStep',t(nt)-t(1));

a = f(:,1);

plot(t, a-Re,'-x')
xlabel('Time (?)')
ylabel('Semi-major axis in [a-Re] (km)')
axis([-inf, inf, -inf, inf])
grid
return

function  dydt = rates(ti,y)
    dfdt = zeros(1,1);
    a = y(1);
        n = sqrt(mu/a^3);
        h = a - Re;
H = 8.4; %km
rho0 = 1.225; %kg/m^3
rho_pe = rho0*exp(-h/H);
    dadt = -BC*rho_pe* 1000*(a^2)*(n/(2*pi))*1*(1+(3/4)*(e^2)*(1));
        fprintf('a\th\n');
        fprintf('%10g %10g \r\n',a,a-Re);

      dydt(1) = dadt;
  end

      function [value,isterminal,direction] = events(t,y)
  % Locate the time when height passes through zero in a decreasing
  % direction
  % and stop integration.
  value = (y(1)-(Re+0));%-(Re+187); % detect a-Re = 0
  isterminal = 1; % stop the integration
  direction = []; % negative direction
      end
  end

Problem solved!!