# Different result with PeriodicCallback

**URL:** <https://discourse.julialang.org/t/different-result-with-periodiccallback/41273>\
**Category:** Modelling & Simulations\
**Tags:** question\
**Created:** [June 12, 2020, 3:46pm UTC](https://discourse.julialang.org/t/different-result-with-periodiccallback/41273 "2020-06-12T15:46:59Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![ohmsweetohm1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ohmsweetohm1/32/49126_2.png) [@ohmsweetohm1](https://discourse.julialang.org/u/ohmsweetohm1)\
**Post date:** [June 12, 2020, 3:46pm UTC](https://discourse.julialang.org/t/different-result-with-periodiccallback/41273/1 "2020-06-12T15:46:59Z")

</div>

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

```julia
using DifferentialEquations, PyPlot

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

f1(du, u, p, t) = du[1] = u.a
f2(du, u, p, t) = du[1] = 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[1] for x in sol1.u], ".-", lw=2, label="solution with callback")
ax1.plot(sol2.t, [x[1] 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.

 ![img](https://global.discourse-cdn.com/julialang/original/3X/7/d/7d07ed388fdcc3de2de8c7bc4e000bca26a21315.png)

What am I doing wrong?

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [June 12, 2020, 4:07pm UTC](https://discourse.julialang.org/t/different-result-with-periodiccallback/41273/2 "2020-06-12T16:07:50Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ohmsweetohm1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ohmsweetohm1/32/49126_2.png) [@ohmsweetohm1](https://discourse.julialang.org/u/ohmsweetohm1)\
**Post date:** [June 12, 2020, 4:25pm UTC](https://discourse.julialang.org/t/different-result-with-periodiccallback/41273/3 "2020-06-12T16:25:52Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [June 12, 2020, 5:51pm UTC](https://discourse.julialang.org/t/different-result-with-periodiccallback/41273/4 "2020-06-12T17:51:06Z")

</div>

> [@ohmsweetohm1](#):
>
> 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.

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [June 12, 2020, 6:10pm UTC](https://discourse.julialang.org/t/different-result-with-periodiccallback/41273/5 "2020-06-12T18:10:58Z")

</div>

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

`f1(du, u, p, t) = du[1] = p[1]`

```julia
function cb!(int) # the callback function
    int.p[1] = 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.

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [June 12, 2020, 6:20pm UTC](https://discourse.julialang.org/t/different-result-with-periodiccallback/41273/6 "2020-06-12T18:20:23Z")

</div>

![image](https://global.discourse-cdn.com/julialang/original/3X/0/d/0daef778fbf7183073c027d258c9d3a7241fa5e7.png)

---

<div class="post-metadata">

**Author:** ![ohmsweetohm1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ohmsweetohm1/32/49126_2.png) [@ohmsweetohm1](https://discourse.julialang.org/u/ohmsweetohm1)\
**Post date:** [June 12, 2020, 6:28pm UTC](https://discourse.julialang.org/t/different-result-with-periodiccallback/41273/7 "2020-06-12T18:28:59Z")

</div>

Interesting! The `DEDataVector` approach is the recommended way, see: [https://docs.sciml.ai/stable/features/diffeq\_arrays/#DEDataArrays-1](https://docs.sciml.ai/stable/features/diffeq_arrays/#DEDataArrays-1)

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [June 12, 2020, 7:56pm UTC](https://discourse.julialang.org/t/different-result-with-periodiccallback/41273/8 "2020-06-12T19:56:12Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![apo383](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/apo383/32/11272_2.png) [@apo383](https://discourse.julialang.org/u/apo383)\
**Post date:** [June 12, 2020, 10:00pm UTC](https://discourse.julialang.org/t/different-result-with-periodiccallback/41273/9 "2020-06-12T22:00:35Z")

</div>

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](https://docs.sciml.ai/v1.8.0/features/callback_functions.html#Example-1:-Bouncing-Ball-1) 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.

---

<div class="post-metadata">

**Author:** ![ohmsweetohm1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ohmsweetohm1/32/49126_2.png) [@ohmsweetohm1](https://discourse.julialang.org/u/ohmsweetohm1)\
**Post date:** [June 12, 2020, 10:58pm UTC](https://discourse.julialang.org/t/different-result-with-periodiccallback/41273/10 "2020-06-12T22:58:18Z")

</div>

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.

---

<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:** [June 13, 2020, 8:24am UTC](https://discourse.julialang.org/t/different-result-with-periodiccallback/41273/11 "2020-06-13T08:24:28Z")

</div>

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:

```julia
using OrdinaryDiffEq, Plots

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

f1(du, u, p, t) = du[1] = u.a
f2(du, u, p, t) = du[1] = 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[1] for x in sol1.u], linestyle = :dash, legend = :topleft, lw=2, label="solution with callback")
plot!(sol2.t, [x[1] for x in sol2.u], linestyle = :dot, lw=5, label="solution without callback")

```

![res](https://global.discourse-cdn.com/julialang/original/3X/b/d/bdbdbc4790b526851ffb280ad8b19f51d1b5646f.png)

---

<div class="post-metadata">

**Author:** ![ohmsweetohm1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ohmsweetohm1/32/49126_2.png) [@ohmsweetohm1](https://discourse.julialang.org/u/ohmsweetohm1)\
**Post date:** [June 13, 2020, 9:24am UTC](https://discourse.julialang.org/t/different-result-with-periodiccallback/41273/12 "2020-06-13T09:24:03Z")

</div>

Thank you!

---

<div class="post-metadata">

**Author:** ![ohmsweetohm1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ohmsweetohm1/32/49126_2.png) [@ohmsweetohm1](https://discourse.julialang.org/u/ohmsweetohm1)\
**Post date:** [June 14, 2020, 4:35pm UTC](https://discourse.julialang.org/t/different-result-with-periodiccallback/41273/13 "2020-06-14T16:35:54Z")

</div>

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

```julia
using DifferentialEquations, PyPlot, Sundials

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

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

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[1] = 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[1] for x in sol1.u], ".-", lw=2, label="solution with callback")
ax1.plot(sol2.t, [x[1] for x in sol2.u], ".-", lw=2, label="solution without callback")
ax1.plot(sol1.t, trueSol.(sol1.t), lw=2, "--", label="true solution")
legend()

```

 ![img](https://global.discourse-cdn.com/julialang/original/3X/b/c/bcb60bd9e6a3a5f65aa8fec19985b3cbb755a38d.png)

---

<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:** [June 14, 2020, 4:47pm UTC](https://discourse.julialang.org/t/different-result-with-periodiccallback/41273/14 "2020-06-14T16:47:34Z")

</div>

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.