Wrong derivatives saved with SavingCallback

General Usage
, ,
1

I’m trying to save the derivative values from an ODE using SavingCallback from DiffEqCallbacks.jl. Why is manual_du not the same as save_du.saveval in this code? Am I misunderstanding when these values are saved, or is this is an issue?

using OrdinaryDiffEq, DiffEqCallbacks
function f!(du, u, p, t)
    du[1] = 0.3u[1]
    du[2] = 0.2u[1] + u[2]
    return nothing
end
save_function = (u, t, integrator) -> (_u = similar(u); get_du!(_u, integrator); return _u)
save_du = SavedValues(Float64, Vector{Float64})
t0, t1 = 0.0, 1.0
saveat = t0:0.1:t1
cb = SavingCallback(save_function, save_du; saveat=saveat)
prob = ODEProblem(f!, [1.0,1.0], (t0,t1), callback=cb)
sol = solve(prob, Tsit5(), saveat=saveat)

manual_du = [zeros(2) for _ in eachindex(saveat)]
for j in eachindex(sol)
    u = sol.u[j] 
    f!(manual_du[j], u, nothing, nothing)
end 
manual_du[1] - [0.3 * 1.0, 0.2 * 1.0 + 1.0]
save_du.saveval[1] - [0.3 * 1.0, 0.2 * 1.0 + 1.0]
manual_du - save_du.saveval

julia> manual_du[1] - [0.3 * 1.0, 0.2 * 1.0 + 1.0]
2-element Vector{Float64}:
 0.0
 0.0

julia> save_du.saveval[1] - [0.3 * 1.0, 0.2 * 1.0 + 1.0]
2-element Vector{Float64}:
 0.009350154485183504
 0.13582658198324649

julia> manual_du - save_du.saveval
11-element Vector{Vector{Float64}}:
 [-0.009350154485183504, -0.13582658198324649]
 [-0.0002137942991308761, -0.0032172187877890224]
 [-0.015279228394231514, -0.26106071340989256]
 [-0.005577907363575996, -0.09867494524327491]
 [-0.0320137675756863, -0.6656221051911233]
 [-0.021712550094837635, -0.46684142492757896]
 [-0.01109761449724933, -0.2469357592801935]
 [-0.00015940524436997006, -0.0036734640814164443]
 [-0.023582896950197274, -0.6267856404645369]
 [-0.011968307129786249, -0.32916801361346115]
 [0.0, 0.0]
2

It’s pulled from the interpolation so it’s only to tolerance. Lower the tolerance on the solve and it should converge.

1 Like