# SavingCallback influences result of ODE problem

**URL:** <https://discourse.julialang.org/t/savingcallback-influences-result-of-ode-problem/91731>\
**Category:** Modelling & Simulations\
**Tags:** diffeqcallbacks\
**Created:** [December 16, 2022, 12:20pm UTC](https://discourse.julialang.org/t/savingcallback-influences-result-of-ode-problem/91731 "2022-12-16T12:20:05Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![Pitaevskii](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pitaevskii/32/35020_2.png) [@Pitaevskii](https://discourse.julialang.org/u/Pitaevskii)\
**Post date:** [December 16, 2022, 12:20pm UTC](https://discourse.julialang.org/t/savingcallback-influences-result-of-ode-problem/91731/1 "2022-12-16T12:20:05Z")

</div>

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.

```julia
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[1]
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[1], tObs[end]), param, callback=cb_save) # Define ODE problem & solve
sol=solve(prob,Vern8(),dt=dt0,adaptive=true,saveat=tObs)

```

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.  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/7/57edec9c35fef512bf15fa1a964e3bbfd36ba01a.png)

I find it counter-intuitive that simply storing a parameter changes the solving algorithm. Am I missing something trivial here?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 16, 2022, 1:45pm UTC](https://discourse.julialang.org/t/savingcallback-influences-result-of-ode-problem/91731/2 "2022-12-16T13:45:53Z")

</div>

> [@Pitaevskii](#):
>
> I find it counter-intuitive that simply storing a parameter changes the solving algorithm. Am I missing something trivial here?

> **[ODE Solvers · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Explicit-Runge-Kutta-Methods)**
>
> Documentation for DifferentialEquations.jl.

> 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`.

---

<div class="post-metadata">

**Author:** ![Pitaevskii](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pitaevskii/32/35020_2.png) [@Pitaevskii](https://discourse.julialang.org/u/Pitaevskii)\
**Post date:** [December 16, 2022, 2:52pm UTC](https://discourse.julialang.org/t/savingcallback-influences-result-of-ode-problem/91731/3 "2022-12-16T14:52:02Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 16, 2022, 2:57pm UTC](https://discourse.julialang.org/t/savingcallback-influences-result-of-ode-problem/91731/4 "2022-12-16T14:57:42Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Pitaevskii](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pitaevskii/32/35020_2.png) [@Pitaevskii](https://discourse.julialang.org/u/Pitaevskii)\
**Post date:** [December 16, 2022, 3:28pm UTC](https://discourse.julialang.org/t/savingcallback-influences-result-of-ode-problem/91731/5 "2022-12-16T15:28:26Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 16, 2022, 5:31pm UTC](https://discourse.julialang.org/t/savingcallback-influences-result-of-ode-problem/91731/6 "2022-12-16T17:31:34Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Pitaevskii](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pitaevskii/32/35020_2.png) [@Pitaevskii](https://discourse.julialang.org/u/Pitaevskii)\
**Post date:** [December 20, 2022, 2:31pm UTC](https://discourse.julialang.org/t/savingcallback-influences-result-of-ode-problem/91731/7 "2022-12-20T14:31:13Z")

</div>

In this case `p` remains unchanged, I simply store the value `p[1]` 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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 22, 2022, 4:11pm UTC](https://discourse.julialang.org/t/savingcallback-influences-result-of-ode-problem/91731/8 "2022-12-22T16:11:07Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 24, 2022, 6:33am UTC](https://discourse.julialang.org/t/savingcallback-influences-result-of-ode-problem/91731/9 "2022-12-24T06:33:02Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Pitaevskii](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pitaevskii/32/35020_2.png) [@Pitaevskii](https://discourse.julialang.org/u/Pitaevskii)\
**Post date:** [January 3, 2023, 1:27pm UTC](https://discourse.julialang.org/t/savingcallback-influences-result-of-ode-problem/91731/10 "2023-01-03T13:27:30Z")

</div>

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:

```julia
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[1]=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[1]+= -0.5γ*u[1]
end

function affect_ic!(integrator) # Callback functions
    integrator.p[1]+=integrator.p[end]
end
save_func(u,t,integrator)=integrator.p[1]
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[1], tObs[end]), param, callback=cb)
sol=solve(prob,Vern8(),reltol=1e-8,abstol=1e-8,adaptive=true,saveat=tObs)

using Plots
plot(tObs,abs2.(sol[1,:])/n0)

```

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.