Dense solution is a 1st order linear interpolation, whereas 4th order spline was promised

[samuela]

In the docs, it says that Tsit5 will return a "free 4th order interpolant". The expectation is that this spline will have knot points at every intermediate step in the RK solve, etc. However, that's not what happens in reality:

import DifferentialEquations: Tsit5, ODEProblem, solve

T = 1.0
sol = solve(
    ODEProblem((x, _, _) -> x, [1.0], (0, T)),
    Tsit5(),
    u0 = [1.0],
    saveat = (0, T),
)

Then,

julia> include("interp_bug.jl")
retcode: Success
Interpolation: 1st order linear
t: 2-element Array{Float64,1}:
 0.0
 1.0
u: 2-element Array{Array{Float64,1},1}:
 [1.0]
 [2.718281708773344]

And in fact we can see that the interpolation is linear, and only between the beginning and end points:

julia> sol.(0:0.1:T)
11-element Array{Array{Float64,1},1}:
 [1.0]
 [1.1718281708773344]
 [1.343656341754669]
 [1.5154845126320031]
 [1.6873126835093377]
 [1.859140854386672]
 [2.0309690252640067]
 [2.202797196141341]
 [2.3746253670186754]
 [2.5464535378960096]
 [2.718281708773344]

julia> exp.(0:0.1:T)
11-element Array{Float64,1}:
 1.0
 1.1051709180756477
 1.2214027581601699
 1.3498588075760032
 1.4918246976412703
 1.6487212707001282
 1.8221188003905089
 2.0137527074704766
 2.225540928492468
 2.45960311115695
 2.718281828459045

Also, since the interpolant is used by default in the adjoint solve, does this imply that gradients are incorrect as well?

[samuela]

Update: curiously removing the saveat argument returns a proper interpolation.

[ChrisRackauckas]

What you're seeing is documented. However, I realized that the documentation was poorly written here, so I gave it a refresh:

https://docs.sciml.ai/latest/basics/common_solver_opts/

* `dense`: Denotes whether to save the extra pieces required for dense (continuous)
  output. Default is `save_everystep && !isempty(saveat)` for algorithms which have 
  the ability to produce dense output, i.e. by default it's `true` unless the user
  has turned off saving on steps or has chosen a `saveat` value. If `dense=false`, 
  the solution still acts like a function, and `sol(t)` is a linear interpolation 
  between the saved time points.
* `saveat`: Denotes specific times to save the solution at, during the solving
  phase. The solver will save at each of the timepoints in this array in the
  most efficient manner available to the solver. If only `saveat` is given, then
  the arguments `save_everystep` and `dense` are `false` by default.   
  If `saveat` is given a number, then it will automatically expand to
  `tspan[1]:saveat:tspan[2]`. For methods where interpolation is not possible,
  `saveat` may be equivalent to `tstops`. The default value is `[]`.

Hopefully that explains it? The solution tries to warn you by saying:

Interpolation: 1st order linear

Also, since the interpolant is used by default in the adjoint solve, does this imply that gradients are incorrect as well?

No. It really depends on the adjoint. When the adjoint is being calculated, it knows that QuadratureAdjoint and InterpolatingAdjoint (when checkpointing is false) needs the interpolation, and so it removes saveat/save_everystep to build the interpolant and then uses that for the backsolve and then corrects the result. When InterpolatingAdjoint is used with checkpointing, the dense solution is not needed so it doesn't perform these corrections. With Backsolve it's not needed so it's like checkpointing. With reverse mode adjoints (TrackerAdjoint, ReverseDiffAdjoint, ZygoteAdjoint) it's not changed either since that's not required for discretize-than-optimize approaches. So it's a bit more subtle but it is handled correctly.

Let me know if you need anything else. Or if you have suggestions to make the docs more clear here.

[samuela]

Ah, I see now the documentation that lays this out. Still, I found the behavior quite surprising. I guess the semantics I was expecting would've been something like: No matter what the integrator always needs to build some kind of interpolation/internal representation. If you additionally want to calculate the values at saveat points, it'll do that for you too. But the saveat functionality and interpolation would be represented and stored separately.

Happy to hear that this doesn't affect the adjoints either way though!

[ChrisRackauckas]

We had it that way at one point, but it confused most people. 99% of the time, if you use saveat then you just added the values you wanted to save at and therefore you don't need dense around. But, a lot of people were then having issues with running out of memory, and the reason was because they weren't expecting dense around. Dense is by far the most memory intensive thing you can do, so if you don't need it then it should go away. For this reason, all of our tutorials used to do things like saveat=0.1,save_everystep=false,dense=false, at which point we realized that if you do saveat then you really just want to save those points, at least most of the time.

[JianghuiDu]

So this interpolation is only used in the final output to calculate the results at saveat, and not used in the internal ODE solver, right? My CVOD_BDF code used to return 3rd order Hermite but recently the same code only returns 1st order linear. Is the interpolation method a parameter that can be set in the ODE problem? Can't find information in the doc.

[ChrisRackauckas]

The rules have been the same at least since 2017 IIRC. If you have a case that doesn't follow the documented rules, please open an issue. If you set saveat, the interpolation will be turned off and you'll have linear interpolation.

So this interpolation is only used in the final output to calculate the results at saveat, and not used in the internal ODE solver, right?

And it's used for sol(t). It's used internally for ContinuousCallback event handling, and sometimes used for stage predictors in stiff ODE solvers.