Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Solve seems to miss timesteps specified by PresetTimeCallbacks #805

Closed
winkmal opened this issue Oct 4, 2021 · 4 comments
Closed

Solve seems to miss timesteps specified by PresetTimeCallbacks #805

winkmal opened this issue Oct 4, 2021 · 4 comments

Comments

@winkmal
Copy link

winkmal commented Oct 4, 2021

I am trying to implement feeding pulses in the right-hand side of an ODE system. Let's say times and values are stored in a CSV file called Inputarray_Odm2prod_2h_3d.csv:

3300,   430, 32.64,        32.64
3521,   0,   32.64,        32.64
3660,   60,  249.83539555, 58.001892497
3828,   0,   249.83539555, 58.001892497
10860,  60,  249.83539555, 58.001892497
11028,  0,   249.83539555, 58.001892497
18060,  60,  249.83539555, 58.001892497
18228,  0,   249.83539555, 58.001892497
25260,  60,  249.83539555, 58.001892497
25428,  0,   249.83539555, 58.001892497
32460,  60,  249.83539555, 58.001892497
32628,  0,   249.83539555, 58.001892497
39660,  60,  249.83539555, 58.001892497
39828,  0,   249.83539555, 58.001892497
46860,  60,  249.83539555, 58.001892497
47028,  0,   249.83539555, 58.001892497
54060,  60,  249.83539555, 58.001892497
54228,  0,   249.83539555, 58.001892497
61260,  60,  249.83539555, 58.001892497
61428,  0,   249.83539555, 58.001892497
68460,  60,  249.83539555, 58.001892497
68628,  0,   249.83539555, 58.001892497
75660,  60,  249.83539555, 58.001892497
75828,  0,   249.83539555, 58.001892497
82860,  60,  249.83539555, 58.001892497
83028,  0,   249.83539555, 58.001892497
89700,  430, 32.64,        32.64
89921,  0,   32.64,        32.64
90060,  60,  249.83539555, 58.001892497
90228,  0,   249.83539555, 58.001892497
97260,  60,  249.83539555, 58.001892497
97428,  0,   249.83539555, 58.001892497
104460, 60,  249.83539555, 58.001892497
104628, 0,   249.83539555, 58.001892497
111660, 60,  249.83539555, 58.001892497
111828, 0,   249.83539555, 58.001892497
118860, 60,  249.83539555, 58.001892497
119028, 0,   249.83539555, 58.001892497
126060, 60,  249.83539555, 58.001892497
126228, 0,   249.83539555, 58.001892497
133260, 60,  249.83539555, 58.001892497
133428, 0,   249.83539555, 58.001892497
140460, 60,  249.83539555, 58.001892497
140628, 0,   249.83539555, 58.001892497
147660, 60,  249.83539555, 58.001892497
147828, 0,   249.83539555, 58.001892497
154860, 60,  249.83539555, 58.001892497
155028, 0,   249.83539555, 58.001892497
162060, 60,  249.83539555, 58.001892497
162228, 0,   249.83539555, 58.001892497
169260, 60,  249.83539555, 58.001892497
169428, 0,   249.83539555, 58.001892497
176100, 430, 32.64,        32.64
176321, 0,   32.64,        32.64
176460, 60,  249.83539555, 58.001892497
176628, 0,   249.83539555, 58.001892497
183660, 60,  249.83539555, 58.001892497
183828, 0,   249.83539555, 58.001892497
190860, 60,  249.83539555, 58.001892497
191028, 0,   249.83539555, 58.001892497
198060, 60,  249.83539555, 58.001892497
198228, 0,   249.83539555, 58.001892497
205260, 60,  249.83539555, 58.001892497
205428, 0,   249.83539555, 58.001892497
212460, 60,  249.83539555, 58.001892497
212628, 0,   249.83539555, 58.001892497
219660, 60,  249.83539555, 58.001892497
219828, 0,   249.83539555, 58.001892497
226860, 60,  249.83539555, 58.001892497
227028, 0,   249.83539555, 58.001892497
234060, 60,  249.83539555, 58.001892497
234228, 0,   249.83539555, 58.001892497
241260, 60,  249.83539555, 58.001892497
241428, 0,   249.83539555, 58.001892497
248460, 60,  249.83539555, 58.001892497
248628, 0,   249.83539555, 58.001892497
255660, 60,  249.83539555, 58.001892497
255828, 0,   249.83539555, 58.001892497

The following code works, but is quite inelegant and inefficient:

using DifferentialEquations
using DelimitedFiles
using Plots

function odm2prod(dx, x, adjParams, t, fixedParams)
	V_liq = fixedParams
	k_1, f_1 = adjParams

	rho_1 = k_1*x[1]
    # Differential Equations
	dx[1] = q_in/V_liq*(X_in - x[1]) - rho_1
	dx[2] = q_in/V_liq*(Y_in - x[2])
end

global q_in = 0
global X_in = 32.64
global Y_in = 32.64
x0          = [11.5, 11.55*0.2]
tspan 	= (0.0, 4.0)
const fixedParams = 159
adjParams = [0.19, 29.0, 0]
prob = ODEProblem((dx, x, adjParams, t) -> odm2prod(dx, x, adjParams,
	t, fixedParams), x0, tspan, adjParams)

inputArray = readdlm("Inputarray_Odm2prod_2h_3d.csv", ',', Float64)
dosetimes  = inputArray[:,1]/86400

function affect!(integrator)
	ind_t = findall(t -> t==integrator.t, dosetimes)
	global q_in = inputArray[ind_t[1], 2]
	global X_in = inputArray[ind_t[1], 3]
	global Y_in = inputArray[ind_t[1], 4]
end
cb = PresetTimeCallback(dosetimes, affect!, save_positions=(false, false))
sol = solve(prob, Tsit5(), callback=cb, saveat=1/24)

plot(sol, vars=[1, 2])

with the resulting plot:
grafik

Since I wanted to take advantage of some advanced features of ModelingToolkit.jl, I recoded the model into this version (using Interpolations.jl), which is supposed to do the same, but more efficient:

using ModelingToolkit, DifferentialEquations
using DelimitedFiles
using Plots
using Interpolations

# %% Define RHS of diff. equation
@parameters k_1, f_1, V_liq
@variables t, rate(t), x(t), q_out(t)  
@variables q_in(t), X_in(t), Y_in(t)
D   = Differential(t)           # define an operator for the differentiation w.r.t. time

# %% Read input from CSV file
inputArray = readdlm("Inputarray_Odm2prod_2h_3d.csv", ',', Float64)
inputArray[:,1] /= 86400                                 # Seconds to days
# does nearest-neighbor interpolation! Only 2D
itp_q_in(t)    = extrapolate(interpolate((inputArray[:,1],), inputArray[:,2], Gridded(Constant{Previous}())), 0.0)(t)
itp_X_in(t)    = extrapolate(interpolate((inputArray[:,1],), inputArray[:,3], Gridded(Constant{Previous}())), 0.0)(t)
itp_Y_in(t)    = extrapolate(interpolate((inputArray[:,1],), inputArray[:,4], Gridded(Constant{Previous}())), 0.0)(t)

@register itp_q_in(t) 
@register itp_X_in(t) 
@register itp_Y_in(t) 

@named odm2prod = ODESystem([
    # Rate equations
    rate[1]  ~ k_1*x,
    # Differential Equations
    D(x    ) ~ q_in/V_liq*(X_in - x    ) - rate[1],
    D(q_out) ~ q_in/V_liq*(Y_in - q_out), 
    D(q_in ) ~ 0,  
    D(X_in ) ~ 0, 
    D(Y_in ) ~ 0] 
)

# %% Define initial conditions, parameters
x0          = [11.5, 11.55*0.2, 0, 0, 0]
tspan       = (0.0, 4.0)
const fixedParams = 159
adjParams = [0.19, 29.0]
# Construct ODEProblem
prob    = ODEProblem(structural_simplify(odm2prod), 
            [
                x       => x0[1],
                q_out   => x0[2],
                q_in    => x0[3],       
                X_in    => x0[4],
                Y_in    => x0[5]
            ],
            tspan,
            [
               k_1      => adjParams[1],
               f_1      => adjParams[2],
               V_liq    => fixedParams[1] 
            ]
        )

# %% Construct callback function
function affect!(integrator)
    integrator.u[3] = itp_q_in(integrator.t)
    integrator.u[4] = itp_X_in(integrator.t)
    integrator.u[5] = itp_Y_in(integrator.t)
end
cb = PresetTimeCallback(inputArray[:,1], affect!, save_positions=(false, false))

# %% Solve ODE, plot results
sol = solve(prob, Tsit5(), callback=cb, saveat=1/24)

plot(sol, vars=[x, q_out])

As one can clearly see, q_in is never reset to 0 between feedings, so it accumulates over time, giving a completely wrong result:
grafik

But evaluating itp_q_in at the specified times interactively gives the correct result! So I am not sure where the error is here. Is it in my code, or some of the involved libraries? Maybe something about tolerances!?

Or maybe there is a better way to implement this whole approach, which avoids this issue!?
@ChrisRackauckas
Copy link
Member

The latter won't be more efficient here. The former is going to be better. For performance, it would more sense to encode those values as parameters and then flip the parameters integrator.p inside of the callback. Either that or making the the former type stable by making those values non-global would be best.

As one can clearly see, q_in is never reset to 0 between feedings, so it accumulates over time, giving a completely wrong result:

I double checked the details here:

function affect!(integrator)
    @show integrator.t, itp_q_in(integrator.t)
    @show integrator.u[3]
    integrator.u[3] = itp_q_in(integrator.t)
    integrator.u[4] = itp_X_in(integrator.t)
    integrator.u[5] = itp_Y_in(integrator.t)
end
cb = PresetTimeCallback(inputArray[:,1], affect!, save_positions=(false, false))

# %% Solve ODE, plot results
sol = solve(prob, Tsit5(), callback=cb)

using Plots
plot(sol)
(0.03819444444444445, 430.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (0.04075231481481482, 430.0)
integrator.u[3] = 430.0
(integrator.t, itp_q_in(integrator.t)) = (0.04236111111111111, 0.0)
integrator.u[3] = 430.0
(integrator.t, itp_q_in(integrator.t)) = (0.044305555555555556, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (0.12569444444444444, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (0.12763888888888889, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (0.20902777777777778, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (0.21097222222222223, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (0.2923611111111111, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (0.29430555555555554, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (0.37569444444444444, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (0.3776388888888889, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (0.45902777777777776, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (0.4609722222222222, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (0.5423611111111111, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (0.5443055555555556, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (0.6256944444444444, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (0.6276388888888889, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (0.7090277777777778, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (0.7109722222222222, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (0.7923611111111111, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (0.7943055555555556, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (0.8756944444444444, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (0.8776388888888889, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (0.9590277777777778, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (0.9609722222222222, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (1.0381944444444444, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (1.0407523148148148, 430.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (1.042361111111111, 0.0)
integrator.u[3] = 430.0
(integrator.t, itp_q_in(integrator.t)) = (1.0443055555555556, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (1.1256944444444446, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (1.1276388888888889, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (1.2090277777777778, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (1.2109722222222221, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (1.292361111111111, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (1.2943055555555556, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (1.3756944444444446, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (1.3776388888888889, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (1.4590277777777778, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (1.4609722222222221, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (1.542361111111111, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (1.5443055555555556, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (1.6256944444444446, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (1.6276388888888889, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (1.7090277777777778, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (1.7109722222222221, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (1.792361111111111, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (1.7943055555555556, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (1.8756944444444446, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (1.8776388888888889, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (1.9590277777777778, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (1.9609722222222221, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (2.0381944444444446, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (2.040752314814815, 430.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (2.042361111111111, 0.0)
integrator.u[3] = 430.0
(integrator.t, itp_q_in(integrator.t)) = (2.0443055555555554, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (2.1256944444444446, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (2.127638888888889, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (2.2090277777777776, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (2.2109722222222223, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (2.292361111111111, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (2.2943055555555554, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (2.3756944444444446, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (2.377638888888889, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (2.4590277777777776, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (2.4609722222222223, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (2.542361111111111, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (2.5443055555555554, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (2.6256944444444446, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (2.627638888888889, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (2.7090277777777776, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (2.7109722222222223, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (2.792361111111111, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (2.7943055555555554, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (2.8756944444444446, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (2.877638888888889, 60.0)
integrator.u[3] = 0.0
(integrator.t, itp_q_in(integrator.t)) = (2.9590277777777776, 0.0)
integrator.u[3] = 60.0
(integrator.t, itp_q_in(integrator.t)) = (2.9609722222222223, 60.0)
integrator.u[3] = 0.0

plot

The values that are supposed to be going to zero are clearly going to zero, and the callback is triggering at all of the times from the CSV. That all checks out as correct to me. The fact that X_in is so huge most of the time is why the derivative is always positive.

@ChrisRackauckas
Copy link
Member 

function affect!(integrator)
	ind_t = findall(t -> t==integrator.t, dosetimes)
	@show integrator.t, inputArray[ind_t[1], 2]
	global q_in = inputArray[ind_t[1], 2]
	global X_in = inputArray[ind_t[1], 3]
	global Y_in = inputArray[ind_t[1], 4]
end
cb = PresetTimeCallback(dosetimes, affect!, save_positions=(false, false))
sol = solve(prob, Tsit5(), callback=cb)
(0.03819444444444445, 430.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.04075231481481482, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.04236111111111111, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.044305555555555556, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.12569444444444444, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.12763888888888889, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.20902777777777778, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.21097222222222223, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.2923611111111111, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.29430555555555554, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.37569444444444444, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.3776388888888889, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.45902777777777776, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.4609722222222222, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.5423611111111111, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.5443055555555556, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.6256944444444444, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.6276388888888889, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.7090277777777778, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.7109722222222222, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.7923611111111111, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.7943055555555556, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.8756944444444444, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.8776388888888889, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.9590277777777778, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (0.9609722222222222, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.0381944444444444, 430.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.0407523148148148, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.042361111111111, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.0443055555555556, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.1256944444444446, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.1276388888888889, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.2090277777777778, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.2109722222222221, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.292361111111111, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.2943055555555556, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.3756944444444446, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.3776388888888889, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.4590277777777778, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.4609722222222221, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.542361111111111, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.5443055555555556, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.6256944444444446, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.6276388888888889, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.7090277777777778, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.7109722222222221, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.792361111111111, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.7943055555555556, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.8756944444444446, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.8776388888888889, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.9590277777777778, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (1.9609722222222221, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.0381944444444446, 430.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.040752314814815, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.042361111111111, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.0443055555555554, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.1256944444444446, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.127638888888889, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.2090277777777776, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.2109722222222223, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.292361111111111, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.2943055555555554, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.3756944444444446, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.377638888888889, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.4590277777777776, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.4609722222222223, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.542361111111111, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.5443055555555554, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.6256944444444446, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.627638888888889, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.7090277777777776, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.7109722222222223, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.792361111111111, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.7943055555555554, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.8756944444444446, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.877638888888889, 0.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.9590277777777776, 60.0)
(integrator.t, inputArray[ind_t[1], 2]) = (2.9609722222222223, 0.0)

You can see that the two dosing schemes are not the same. The interpolation gives: (integrator.t, itp_q_in(integrator.t)) = (0.04075231481481482, 430.0), while the direct indexing gives: (integrator.t, inputArray[ind_t[1], 2]) = (0.04075231481481482, 0). So it looks like you didn't use the interpolation library correctly. I haven't used Interpolations.jl so I don't know why it would do nearest in that way. My guess is that it's not nearest but something like <, so being off by floating point error amount can make it hit the wrong value. Whatever is going on, that's an interpolation issue, not a DiffEq issue.

@winkmal
Copy link
Author

winkmal commented Nov 30, 2021

Dear Chris, thanks for looking into this.
Til now, I did not know the @show macro, seems very useful for debugging.

I found that when I add a very small number to dosetimes vector in the PresetTimeCallback call, it actually works, like:
cb = PresetTimeCallback(inputArray[:,1].+1e-6, affect!, save_positions=(false, false))
grafik

The reason I went with Interpolations.jl in the first place was actually their Gridded(Constant{Previous}() mode.
Maybe DataInterpolations.jl can do the same, but with better performance/without the roundoff error!?

@ChrisRackauckas
Copy link
Member

DataInterpolations.jl won't have the round-off issue as it will be exact at the data points. As for performance, I don't know.

@winkmal winkmal mentioned this issue Dec 15, 2021
Open
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@ChrisRackauckas @winkmal