Reputation: 920
I have an ODE to be integrated with ode113
. Inside the function that describes the derivatives I compute some variables that I would like to "save" (either store into a file or into an array). How can I do that?
Here is an example (adapted from the documentation). Run the following script:
global number_of_calls
number_of_calls = 0;
[t, y] = ode113(@myfunc, [0 10], 0);
plot(t,y,'-o')
fprintf('Number of calls of ''myfunc'': %d\n', number_of_calls)
fprintf('Number of elements of ''t'' : %d\n', length(t))
which calls the function containing the ODE to be integrated:
function dydt = myfunc(t, y)
global number_of_calls
dydt = 2*t;
a = 2*y; % <--- how can I save a?
number_of_calls = number_of_calls + 1;
end
I thought of writing to file the values of a
from inside myfunc
, but as you can see myfunc
is called multiple times before actually storing the current value for the independent variable, t
, and in the end I would have two arrays a
and t
of different size.
Upvotes: 1
Views: 159
Reputation: 125854
Here's one way you can write myfunc
to do this. It saves the results to a global vector a
, which must be set to []
before each call to ode113
. It uses a persistent variable to keep track of the previous time step:
function dydt = myfunc(t, y)
dydt = 2*t;
temp = 2*y; % Extra value to save
global a
persistent tPrevious
if isempty(a) % Initialize
a = temp;
elseif (t > tPrevious) % Moving to a new time step, so add an element
a = [a; temp];
else % Refining time step downwards or recomputing at the same time step
a(end) = temp;
end
tPrevious = t; % Store prior time step
end
And you can check that it's working like so:
global a;
a = [];
[t, y] = ode113(@myfunc, [0 10], 0);
isequal(2.*y, a)
ans =
logical
1 % a exactly equals 2.*y
Upvotes: 1