Callback problem with DifferentialEquations

I am stuck with a problem considering a SavingCallback from the DifferentialEquations suite.

I could not create a MWE, thus I will explain the problem.
I have an ODEProblem where I use a SavingCallback for the calculation of some internals. In the calculation I need the du vector, thus I have created a callback function (I have seen that on stackoverflow) similar to

f = (u,t,integrator) -> calculate_internals(get_du(integrator), u, t)

If I solve this, I get an UndefRefError: access to undefined reference error, because it tries to access the fsallast property that apparently does not exist. However, if I check with propertynames(integrator) it does appear in the list (I checked this in the Juno debugger). Is this a bug?

Note that I can solve the problem without providing the callback or if I instead of providing the du vector in the callback function, I provide a random vector with the same size.

Are you calling this in initialize? I need a bit more information than this. Share some code.

Thanks for getting back to me!

What do you mean with initialize?

Originally, the calculate_internals function is called within the ode function. If I want to extract the internals, I use the calculate_internals function in a SavingCallback. It worked when I did not use the get_du function previously and as said it works if I supply a random vector with the length of the states.

I don’t have access to the code atm.

I was thinking about augmenting the states and have all the internals as additional states by setting du=0 and setting the states u within the ode function. Is there any speed compromise that comes with it, or any other disadvantages?

Yeah I think the issue is just when it’s called at the initial point, If you do save_start = false in the SavingCallback it should be fine though.

1 Like

save_start = false does not work.

I am using now integrator(t, Val{1}) as you suggested in https://github.com/SciML/OrdinaryDiffEq.jl/issues/785#issuecomment-504962062. That works.

Unfortunately, the approach does not give the same derivative values as the get_du approach.

However, I have managed to create a MWE (partly taken from the official docs):

using DifferentialEquations

function lorenz!(du,u,p,t)
 du[1] = 10.0*(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end

u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan)

integ1 = init(prob, Tsit5())
get_du(integ1) # works
integ2 = init(prob, Rodas4P())
get_du(integ2) # throes the mentioned error

Seems to work with all the other Rodas implementations apart from the ones with 4th order.

And a MWE in combination with a callback where it fails for all solvers (also taken from the official docs, see here:

using DiffEqCallbacks, OrdinaryDiffEq, LinearAlgebra
prob = ODEProblem((du,u,p,t) -> du .= u, rand(4,4), (0.0,1.0))
saved_values = SavedValues(Float64, Tuple{Float64,Float64})
cb = SavingCallback((u,t,integrator)->(get_du(integrator)[1],norm(u)), saved_values)
sol = solve(prob, Tsit5(), callback=cb) # throws mentioned error

You were right, specifying save_start = false in the SavingCallback did solve the issue if another solver than one from the Rodas4 family is chosen.

I do have another short question:
The time values in the saved_values array do not need to be necessarily in ascending order, is that correct (happens in my case)?

Could you open an issue on OrdinaryDiffEq? This likely comes down to being an issue at the initial point for some customized interpolants.

Done, see here