# SavingCallback influences result of ODE problem

Hello

I have some code for simulating complex differential equations of a certain model that used to work fine, but fails after going from Julia 1.7 to 1.8 and updating all the packages.

``````using DifferentialEquations, DiffEqCallbacks, SparseArrays
function du!(du,u,p,t)                                  # Define differential equations
J,Jm,U,γv,D=p
mul!(D,Jm,u)
@. du = im*J*D - (im*U*abs2(u)+γv)*u
end

save_func(u,t,integrator)=integrator.p
saved_values=SavedValues(Float64,Float64)
cb_save=SavingCallback(save_func,saved_values,saveat=tObs)

M=size(ψ0,1); mat=build(jp);  # Set parameters
param=[jp.J, mat.Jm, jp.U, mat.γv, zeros(ComplexF64,M)]
#---------------------------------------

prob=ODEProblem(du!, ψ0, (tObs, tObs[end]), param, callback=cb_save)   # Define ODE problem & solve
``````

Without the SavingCallback there is no problem and the `Vern8()` solver is more than suffucient with the default tolerances. When the Savingcallback is added, the exact same setup gives incorrect results. In the figure I plot the absolute value of one element of the vector `u` that in this model is supposed to stabilize. I find it counter-intuitive that simply storing a parameter changes the solving algorithm. Am I missing something trivial here?

https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Explicit-Runge-Kutta-Methods

Additionally, the following algorithms have a lazy interpolant:

• `BS5` - Bogacki-Shampine 5/4 Runge-Kutta method. (lazy 5th order interpolant).
• `Vern6` - Verner’s “Most Efficient” 6/5 Runge-Kutta method. (lazy 6th order interpolant).
• `Vern7` - Verner’s “Most Efficient” 7/6 Runge-Kutta method. (lazy 7th order interpolant).
• `Vern8` - Verner’s “Most Efficient” 8/7 Runge-Kutta method. (lazy 8th order interpolant)
• `Vern9` - Verner’s “Most Efficient” 9/8 Runge-Kutta method. (lazy 9th order interpolant)

These methods require a few extra steps in order to compute the high order interpolation, but these steps are only taken when the interpolation is used. These methods when lazy assume that the parameter vector `p` will be unchanged between the moment of the interval solving and the interpolation. If `p` is changed in a ContinuousCallback, or in a DiscreteCallback and the continuous solution is used after the full solution, then set `lazy=false`.

Ok, I didn’t realise this also holds for the SavingCallback. Thanks for the help!

We should probably highlight that then. I’ll make a note of it.

I was just confused because `p` remains unchanged, but I see how it can still affect the solver

I didn’t read your case in detail, I thought you had a `p` modification in there somewhere. If you set `lazy=false` is it fine?

In this case `p` remains unchanged, I simply store the value `p` at points `saveat=tObs`. This is ofcourse not useful, but merely a test to see where the error was. When I put `lazy=false` the simulation runs smoothly again.

Oh I know what’s going on here. Can you open an issue on DiffEqCallbacks so it’s in my queue?

I’m not able to run your code, it’s missing many things in its definitions.

Hi Chris, sorry for the delay. The code I provided indeed doesn’t run, but is just an indicator of where the problem occurred. Here is a working example:

``````using DifferentialEquations, DiffEqCallbacks, SparseArrays, LinearAlgebra

function build(Nn::Int64,BC::String)::Matrix{Float64}
if BC=="open"
Jm=spdiagm(1=>ones(Nn),-1=>ones(Nn)); Jm[1,2]*=2;
elseif BC=="periodic"
Jm=spdiagm(1=>ones(Nn),-1=>ones(Nn)); Jm[1,2]*=2;
Jm[Nn+1,Nn+1]=1.0;
elseif BC=="neumann"
Jm=spdiagm(1=>ones(Nn),-1=>ones(Nn)); Jm[1,2]*=2;
Jm[Nn+1,Nn]*=2;
else
error("No matching boundary condition")
end
return Jm
end
function initialChain(n0::Int64,Nn::Int64,ϕ::Float64,IC::String)::Vector{ComplexF64}
ψ0=zeros(ComplexF64,Nn+1)
if IC=="full" # Central site is full
ψ0[:].=sqrt(n0)*exp(im*ϕ);
elseif IC=="empty" # Central site starts of empty
ψ0[2:Nn+1].=sqrt(n0)*exp(im*ϕ);
ψ0=sqrt(0.1*n0)*exp(im*ϕ);
else
error("Undefined initial condition")
end
return ψ0
end
function du!(du,u,p,t)                  # Differential equations
J,Jm,U,γ,D=p
mul!(D,Jm,u)
@. du = im*J*D - im*U*abs2(u)*u
du+= -0.5γ*u
end

function affect_ic!(integrator)         # Callback functions
integrator.p+=integrator.p[end]
end
save_func(u,t,integrator)=integrator.p
saved_values=SavedValues(Float64,Float64)

U=1; n0=700; Nn=200; J=70; γ=150; ϕ=0.0pi; IC="full"; BC="periodic";
tObs=Array(0:0.001:1); dosetimes=Array(range(0.5/U/40,0.5/U,40))
D=zeros(ComplexF64,Nn+1); Jm=build(Nn,BC); ψ0=initialChain(n0,Nn,ϕ,IC)

if IC=="full"
δJ=(J-100.0)/40
param=[100.0, Jm, U, γ, D, δJ]
elseif IC=="empty"
δJ=(J-40.0)/40
param=[40.0, Jm, U, γ, zeros(ComplexF64,M), δJ]
end
cb_preset=PresetTimeCallback(dosetimes,affect_ic!,save_positions=(false,false))
cb_save=SavingCallback(save_func,saved_values,saveat=tObs)
cb=CallbackSet(cb_preset,cb_save)

prob=ODEProblem(du!, ψ0, (tObs, tObs[end]), param, callback=cb)
Using only the PresetTimeCallback, that changes `p`, works fine. Taking `lazy=false` does not seem to change the result. Using the SavingCallback does require `lazy=false`, even when you leave out the PresetTimeCallback.