Pier Paolo
Pier Paolo

Reputation: 920

How to save a variable from within the integration with ode113?

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

Answers (1)

gnovice
gnovice

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

Related Questions