# Different result with PeriodicCallback

I am trieing to understand the `PeriodicCallback` method and noticed something weird.

``````using DifferentialEquations, PyPlot

mutable struct SimType{T} <: DEDataVector{T}
x::Vector{T}
a::T
end

f1(du, u, p, t) = du = u.a
f2(du, u, p, t) = du = if 10. <= t <= 20.  1. else 0. end

function cb!(int) # the callback function
int.u.a = if 10. <= int.t <= 20.  1. else 0. end
return
end

u0 = SimType([0.], 0.)
tspan = (0., 30.)
prob1 = ODEProblem(f1, u0, tspan)
sol1 = solve(prob1, Tsit5(); callback = PeriodicCallback(cb!, 1., save_positions = (false,false)))

prob2 = ODEProblem(f2, [0.], tspan)
sol2 = solve(prob2, Tsit5())

function trueSol(t)
if t <= 10.
0.
elseif 10. < t <= 20.
-10. + t
else
10.
end
end

pygui(true)
fig = figure(figsize=(12, 5))
ax1 = subplot(111)
ax1.plot(sol1.t, [x for x in sol1.u], ".-", lw=2, label="solution with callback")
ax1.plot(sol2.t, [x for x in sol2.u], ".-", lw=2, label="solution without callback")
ax1.plot(sol1.t, trueSol.(sol1.t), lw=2, "--", label="true solution")
legend()
``````

Here is a picture of the problem.

What am I doing wrong?

I can’t try this right now, but I don’t understand why you have `save_positions = (false,false)`. Have you run it without that?

I’ve never used that option - this is what the manual says:

`save_positions` : Boolean tuple for whether to save before and after the `affect!` . The first save will always occcur (if true), and the second will only occur when an event is detected. For discontinuous changes like a modification to `u` to be handled correctly (without error), one should set `save_positions=(true,true)` ."

Makes no difference for the result. However, if set to `(true, true)` I get double entries in the time vector of the solution.
I also don’t get why a correct computation would depend on this setting.

I don’t either but read that line in the docs and thought it was worth a try. I’m going to try to run your example once some things finish installing and precompiling.

If I make the following modifications (passing the slope in a parameter vector p) I get closer:

`f1(du, u, p, t) = du = p`

``````function cb!(int) # the callback function
int.p = if 10. <= int.t <= 20.  1. else 0. end
return
end
``````

`p = [0.0,]`
`prob1 = ODEProblem(f1, u0, tspan,p)`

Now the slope is correct but it doesn’t stop until 21.0 s. [ change the `<=` to `<` fixes that].

I have not used a custom type like your SimType for the integration variable, but have only used simple vectors.

Interesting! The `DEDataVector` approach is the recommended way, see: https://docs.sciml.ai/stable/features/diffeq_arrays/#DEDataArrays-1

Yes, it would appear you are trying to do it the way the docs say. I don’t know.

I think it’s time to call DifferentialEquationsMan (I mean @ChrisRackauckas).

1 Like

ODE integrators are usually making a polynomial approximation to the true solution. If you have a change in slope at the event, there is no polynomial approximation for that. The `(true, true)` tells the integrator to use different approximations for each side of the event. It no longer feels compelled to make things C1 continuous.

By the way, I would also advise against a discrete event such as `if 10. <= int.t <= 20. 1.`, which gives the integrator no continuous information to triangulate the event. Better to provide a zero-crossing like the bouncing ball demo in docs. A continuous function helps the integrator know when to adjust its steps as it gets closer to the event.

Also by the way, if you do a discrete event, you needn’t do `if ... 1 else 0`, you can just make it logical, `10. <= int.t <= 20.`. But again, better to do a continuous zero crossing.

Thank you for the input but I do not think that this is the issue here. Otherwise why would the version of @klaf using the parameter in the callback work?
In the bouncing ball example the problem is that the event is not detected. When using the periodic callback you already know the time of the event so there should be no problem.

Hey,

1. When using the `DEDataArray`, you have to remember to update all of the cache variables. Notice the docs show the `full_cache` iterator.
2. When solving an ODE with discontinuities, you have to make sure you resolve the discontinuities if you want a high precision answer. Here, you know them in advance, so you can set `tstops = [10.0,20.0]` to hit those points exactly during the integration.

With those two facts together it all works out:

``````using OrdinaryDiffEq, Plots

mutable struct SimType{T} <: DEDataVector{T}
x::Vector{T}
a::T
end

f1(du, u, p, t) = du = u.a
f2(du, u, p, t) = du = if 10. <= t < 20.  1. else 0. end

function cb!(int) # the callback function
for c in full_cache(int)
c.a = if 10.0 <= int.t < 20.  1. else 0. end
end
return
end

u0 = SimType([0.], 0.)
tspan = (0., 30.)
prob1 = ODEProblem(f1, u0, tspan)
sol1 = solve(prob1, Tsit5(); callback = PeriodicCallback(cb!, 1., save_positions = (false,false)))

prob2 = ODEProblem(f2, [0.], tspan)
sol2 = solve(prob2, Tsit5(), tstops = [10.,20.])

function trueSol(t)
if t <= 10.
0.
elseif 10. < t <= 20.
-10. + t
else
10.
end
end

plot(sol1.t, trueSol.(sol1.t), lw=2, label="true solution")
plot!(sol1.t, [x for x in sol1.u], linestyle = :dash, legend = :topleft, lw=2, label="solution with callback")
plot!(sol2.t, [x for x in sol2.u], linestyle = :dot, lw=5, label="solution without callback")
`````` 2 Likes

Thank you!

Could this `DEDataVector`-approach also be applied to DAE problems?
I had no luck so far and went with the parameter-approach.

``````using DifferentialEquations, PyPlot, Sundials

mutable struct SimType{T} <: DEDataVector{T}
x::Vector{T}
a::Real
end

f1(du, u, p, t) = du = u.a
f2(res, du, u, p, t) = res = p - du

function cb1!(int) # the callback function
for c in full_cache(int)
c.a = if 10.0 <= int.t < 20.  1. else 0. end
end
return
end
function cb2!(int) # the callback function
int.p = if 10.0 <= int.t < 20.  1. else 0. end
return
end

u0 = SimType([0.], 0.)
tspan = (0., 30.)
prob1 = ODEProblem(f1, u0, tspan)
sol1 = solve(prob1, Tsit5(); callback = PeriodicCallback(cb1!, 1., save_positions = (false,false)))

u0 = [0.]
du0 = [0.]
p = [0.]
prob2 = DAEProblem(f2, u0, du0, tspan, p, differential_vars = [true])
sol2 = solve(prob2, Sundials.IDA(); callback = PeriodicCallback(cb2!, 1., save_positions = (false,false)))

function trueSol(t)
if t <= 10.
0.
elseif 10. < t <= 20.
-10. + t
else
10.
end
end

pygui(true)
fig = figure(figsize=(12, 5))
ax1 = subplot(111)
ax1.plot(sol1.t, [x for x in sol1.u], ".-", lw=2, label="solution with callback")
ax1.plot(sol2.t, [x for x in sol2.u], ".-", lw=2, label="solution without callback")
ax1.plot(sol1.t, trueSol.(sol1.t), lw=2, "--", label="true solution")
legend()
``````

It definitely can’t be used with Sundials, and generally with implicit it will have undefined behavior. I’m probably removing it from the docs because of these issues, and I wouldn’t quite recommend it anymore.