What are the details of the interpolation scheme done in ode45?

I understand that ode45 uses Dormand-Prince method to generate the solution at the adaptive time steps. So in the case when tspan input is with more than two elements, the documentation points out that it doesn't affect internal time step. So even though it returns the answer at the specified tspan values, it's not guaranteed that the solver stepped exactly to those values but rather somewhere close, since it's adaptive in nature. My question is what does it use to produce the answer at the desired node since the solver didn't step there? What are the details of this interpolation scheme?

This is the documentation found in the help for ode45 (below):

Specifying tspan with more than two elements does not affect the internal time steps that the solver uses to traverse the interval from tspan(1) to tspan(end). (All solvers in the ODE suite obtain output values by means of continuous extensions of the basic formulas) [WHAT DOES THIS SENTENCE MEAN?!]. (Although a solver does not necessarily step precisely to a time point specified in tspan, the solutions produced at the specified time points are of the same order of accuracy as the solutions computed at the internal time points) [WHAT METHOD IS USED FOR THIS? TO GET THE SAME ORDER OF ACCURACY].

Thanks in advance!

Thanks for the help!

Answer 1

[Kelly Kearney]

The interpolation is done by the ntrp45 function (<matlabroot>/toolbox/matlab/funfun/private/ntrp45.m); you can open that and examine.

Answer 2

[Dragan Plakalovic]

Thanks.

I wrote a standard 4-step RK4 method (fixed step), and am finding out that it is slower than ode45. That surprised me. Say I am only integrating for 1 second, 0 to 1, at 0.001 steps...that's 1000 steps of RK4. Now ode45 will do this in only a few internal steps (for the problem I am doing, and the given tolerances), say 10 steps at 0.1 intervals (hypothetically). But if I specified 0:0.001:1, it means that it has to interpolate (using ntrp45) the other 990 points.

This would mean then, I suppose, that the overhead for interpolation, and step refinement based on relative/absolute tolerance is smaller than doing one step of RK4.

Is this plausible to you?

Thanks!

[Kelly Kearney]

Sounds plausible, though it would depend on exactly how you implemented the RK4 function. The Mathworks used to offer fixed-step solvers (ode1, ode2, ode3, ode4, and ode5) on their website through one of the tech notes, but now that they've redone the website I can't find them... anyone know where those are now?

[Dragan Plakalovic]

That'd be nice to see their implementation. Although I think my implementation is fine and very simple. There isn't anything complicated creating a bottleneck. I have a separate function defining the derivative formulation, and calling it in the RK's 4 steps.

Answer 3

[Dragan Plakalovic]

I was able to find Matlab's ode4.m code just through google search and it's just as slow as the code I wrote, if not slower. It basically comes down to slick interpolation technique in the ntrp45.m where efficiency is not sacrificed but you gain on adaptive time stepping. I'll look further into ntrp45 to see if they interpolate every point separately or it's some sort of vectorization involved to make it even faster.

Thanks again for your original reply.

[金灵 周]

Hello, I have the same confusion.whether ntrp45 interpolate every point separately or it's some sort of vectorization involved to make it even faster?May I ask if you are making progress in this question

Thank you for your reply.