Problem with returning to Main function from subfunction

4 次查看（过去 30 天）
David
David 2012-1-11
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
  2 个评论
Andrew Newell
Andrew Newell 2012-1-11
This ODE solver is from MATLAB Central. Have you tried contacting the author? Also, have you considered using the MATLAB solvers like ode45?
David
David 2012-1-11
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.

请先登录，再进行评论。

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

采纳的回答

Walter Roberson
Walter Roberson 2012-1-12
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.

更多回答（2 个）

David
David 2012-1-12
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!
David
David 2012-1-16
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!!

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by