PresetTimeCallback does not reach u[1] if save_positions = false & using interpolations

Specific Domains
,
1

Hey All,

I’m trying to calculate the error of a differential equations sol object with respect to specific timepoints.
The system is made up as follows:

function dAdt(dA, A, p, t)
    CL, V1, Q, V2 = p
    k₁₀ = CL / V1
    k₁₂ = Q / V1
    k₂₁ = Q / V2
    dA[1] = A[2] * k₂₁ - A[1] * (k₁₀ + k₁₂)
    dA[2] = A[1] * k₁₂ - A[2] * k₂₁
end

I am using a PresetTimeCallback to change the system at timepoints which differs between individuals. Essentially I put a dose inside the volume (V1) of compartment 1 (A[1]):

dose_timepoint = 0.
dose = 2000.
affect!(integrator) = integrator.u[1] += (dose / integrator.p[2])
callback = PresetTimeCallback(dose_timepoint, affect!; save_positions=(false, false)

I am using save_positions=(false, false) here since I want the sol object to only return my saveat indexes:

parameters = [200., 2000., 100., 3000.]
measurement_times = [0., 1., 4. 24., 48]
prob = ODEProblem(dAdt, [0., 0.], (-20., 120.), parameters)
sol = solve(prod, tstops=dose_timepoint, saveat=measurement_times, callback=callback)
# for example returns u[1]=[10., 8., 4., 0., 0.]

This works as expected. When I want to use interpolation however the dose from the callback is added to A[1], but this is not reflected in u[1]:

sol = solve(prod, tstops=dose_timepoint, callback=callback)
sol(measurement_times)
# for example returns u[1]=[-5., 3., 1., 0., 0.]

What you can see above is that A[1] is filled with the dose from the callback, causing dA[1] to be negative. This time however u[1] remains 0., causing the result to become negative. After some time the concentration in A[2] is large enough to actually start filling A[1], so A[1] becomes postive again.
Setting save_positions = (true, true) for the callback fixes the issue.

Is this a expected issue when using interpolation? I of course can manage with using saveat, and setting save_positions = true when I need to use interpolation, but I would expect save_positions to only save the positions to the sol.t, and not affect the solution in other ways, as seems to be the case.
What are your thoughts? Is this a bug?

2

Yes, it’s expected. If you don’t save the discontinuity points, how could you ever interpolate correctly? You need to save each of the points and two points at a jump discontinuity in order to do it.

3

Right. I figured it was likely going to be something to do with the way interpolation works.
When using interpolation you don’t really care of course that initially the points are stored in the sol.t output.

Thanks as always for the fast response Chris!

4

This illustration with the bouncing ball should explain it:

using OrdinaryDiffEq, Plots

function f(du,u,p,t)
  du[1] = u[2]
  du[2] = -p
end

function condition(u,t,integrator) # Event when event_f(u,t) == 0
  u[1]
end

function affect!(integrator)
  integrator.u[2] = -integrator.u[2]
end

cb = ContinuousCallback(condition,affect!)

u0 = [50.0,0.0]
tspan = (0.0,15.0)
p = 9.8
prob = ODEProblem(f,u0,tspan,p)
sol = solve(prob,Tsit5(),callback=cb)
plot(sol)
savefig("plot1.png")

cb2 = ContinuousCallback(condition,affect!,save_positions=(false,false))
sol = solve(prob,Tsit5(),saveat=0.75,callback=cb2)
plot(sol)
savefig("plot2.png")

plot1 plot2

Plot 2 isn’t wrong, it’s just interpolating over the point with the jump discontinuity which introduces an error. Plot 1, because it’s saving twice at the discontinuity, is resolving it exactly, giving a nice vertical line. This is why the double save is the default in the callbacks, as it’s the only way to do this precisely. However, if you do want to just save at specific points, you can turn it off and the interpolation will fill it in via a natural way, but of course without exact resolution of the jump you will have some artifacts.

5

Clear explanation, thanks!