SavingCallback influences result of ODE problem

Specific Domains Modelling & Simulations
1

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[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

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

2

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.

3

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

4

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

5

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

6

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?

7

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.

8

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

9

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

10

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