How to make ode45 store/return only results of every n-th time step to reduce memory use?

Hi,

I am trying to run a relatively long simulation with a small time step with ode45, meaning the resulting y matrix is very big. I want to keep the time step small for accuracy, but I do not need the result at each time step. When running a longer simulation the memory becomes a problem so makin the output sparser would be very beneficial.

Thank you!

[Torsten]

@Stephen23

I assumed that with "keeping the time step small" the OP meant that "RelTol" and "AbsTol" in his simulation were chosen small, not that the large output was a consequence of setting the "tspan" vector accordingly. The latter will not enhance precision as you correctly pointed out.

Answer 1

[Torsten]

You can specify the times when the solution should be saved in the input variable "tspan" to "ode45".

And you should test whether "ode15s" can make larger time steps - maybe your system is stiff so that "ode45" is not suited as solver.

Answer 2

[Walter Roberson]

Set your tspan to be a vector of three elements (more correctly, use more than 2 elements). When you call ode45(), set it up to ignore the t and y outputs.

Use odeset to create ode*() options that include OutputFcn . Have the output function use a persistent variable to count how many times it has been invoked. "most" of the time throw away the data. Every N'th time it is called, store the data in an accessible place.

OutputFcn will be called for each successful time step; you ignore most of the time steps and record every N'th.

Meanwhile the internal recording of the outputs will only happen at the times in your tspan, and you will have stripped down your tspan to the minimum needed so that it does not choose to record everything (tspan with 2 entries.)

[Torsten]

@Walter Roberson

Seems OutputFcn is only called at the three times specified in "tspan". Or is there something wrong with the example code ?

f = @(t,y) exp(-0.05*t)*sin(t);
tspan = [0 0.5 1];
options = odeset('OutputFcn',@myOutputFcn);
y0 = 0;
[T,Y] = ode45(f,tspan,y0,options);

t = 1×2

0 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

t = 0.5000

t = 1

t =

     []

function status = myOutputFcn(t,y,flag)
  status = 0;
  t
end

[Mike Hosea]

Clearly the interface is designed to require the user to supply the time points in situations like this, but there's no adaptivity of the output points in that case.

I don't know if the need for this is common enough to merit supporting it more direclty, but I could imagine an enhancement to support it directly, albeit in a slightly kludgy way, if refine < 0 is the number of consecutive steps to throw away.

[Torsten]

In this case, I don't understand what "OutputFcn" is good for if the function is only called at the times specified in "tspan" (assuming "tspan" is a vector with more than two elements). Postprocessing the results of the integration could be done easily after the simulation has finished.

[Mike Hosea]

Torsten is right. I was thinking that outputFcn was called after every successful step. In fact, it is called after every successful step in which there are any output points. If there are none, it will not get called.

If the problem isn't stiff, I guess you could minimize the number of time steps by using the highest order RK method available, i.e. ode89 and setting Refine=1.

Other than that, I think you have to specify output points up front to minimize storage.

[Torsten]

@Stephen23

What time points do you imagine that OUTPUTFCN shoudl be called at: all function evaluations including the invalid solutions and backtracking?

"OutputFcn" should be called after each successful integration step - independent of what is set in "tspan".

[Mike Hosea]

@Stephen23, I was assuming that the user was using length(tspan) == 2 and then complaining that the number of output points was too large. I thought their reference to the time step was made with respect to tolerance setting. As I read your comment, I see the other interpretation that you suspected. It is plausible, but I don't know which interpretation is correct at this point.

[Stephen23]

@Torsten: isn't your suggested behavior already provided by specifying TSPAN with only the end points?

[Torsten]

Isn't your suggested behavior already provided by specifying TSPAN with only the end points?

No, it doesn't solve the OP's problem. Specifying TSPAN with only the end points means that all the solutions at successful integration steps will be written in the Y-matrix returned by the ode-solver - and this matrix can become very large.

In my opinion, plotting can be done after the solution is returned by the ODE-solver.

[Walter Roberson]

"OutputFcn" should be called after each successful integration step

That's how it is documented.

- independent of what is set in "tspan".

Unfortunately the documentation does not specify the relationship between recording output and OutputFcn calls. I think it is a valid interpretation to assume that OutputFcn is called after every successful integration step regardless of whether the step is one to be recorded in the output arrays. That might turn out to be inaccurate in practice but I think it is a valid reading.

[Torsten]

For the given example, OutputFcn is called for the three times specified in "tspan". You can see the difference if you replace "tspan" by the vector that only contains the endpoints of the integration interval.

f = @(t,y) exp(-0.05*t)*sin(t);
options = odeset('OutputFcn',@myOutputFcn);
y0 = 0;

tspan = [0 0.5 1];
[T,Y] = ode45(f,tspan,y0,options);

t = 1×2

0 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

t = 0.5000

t = 1

t =

     []

tspan = [0 1];
[T,Y] = ode45(f,tspan,y0,options);

t = 1×2

0 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

t = 1×4

0.0250 0.0500 0.0750 0.1000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

●

t = 1×4

0.1250 0.1500 0.1750 0.2000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

t = 1×4

0.2250 0.2500 0.2750 0.3000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

t = 1×4

0.3250 0.3500 0.3750 0.4000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

t = 1×4

0.4250 0.4500 0.4750 0.5000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

t = 1×4

0.5250 0.5500 0.5750 0.6000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

t = 1×4

0.6250 0.6500 0.6750 0.7000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

t = 1×4

0.7250 0.7500 0.7750 0.8000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

t = 1×4

0.8250 0.8500 0.8750 0.9000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

t = 1×4

0.9250 0.9500 0.9750 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

t =

     []

T

T = 41×1

0 0.0250 0.0500 0.0750 0.1000 0.1250 0.1500 0.1750 0.2000 0.2250

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function status = myOutputFcn(t,y,flag)
  status = 0;
  t
end

[Mike Hosea]

Yes, it turns out that the outputFcn is only called if there is a non-empty set of outputs to transmit to it. I guess that is in keeping with it being literally an "output function" and not, more generically, a step callback function. Probably most use cases involve length(tspan)==2.

I have reported the documentation issue.

Answer 3

[Mike Croucher]

I would suggest trying the new ode interface The new solution framework for Ordinary Differential Equations (ODEs) in MATLAB R2023b » The MATLAB Blog - MATLAB & Simulink

One approach that might work for you is the new solutionFcn. This is covered in the blog post above but, for completeness, here's the example I used

myODE = ode(ODEFcn=@(t,y) exp(-0.05*t)*sin(t),InitialTime=0,InitialValue=0); % Define the ODE problem
solFcn = solutionFcn(myODE,0,20)                                            % Use it to create a solutionFcn                                

solFcn = function_handle with value:
    @ode.solutionFcn/interpolate

I can now evaluate this wherever I like within the interval defined in the solutionFcn. E.g.

solFcn(15)

ans = 1.3404

This way you can keep only the exact time points you want.

[Torsten]

@Mike Croucher

The OP's problem is that the ode integrator internally saves the solution for too many time steps such that available memory becomes low during a computation. I think this problem is not solved by your suggestion. Somewhere, the solution of the ODE must be saved at discrete times in order to make interpolation using the function handle "solFcn" after the computation has finished - and I don't see that this solution is somehow reduced in size by the code.

[Mike Croucher]

Yes, you are right! Sorry for the noise