Different result with PeriodicCallback

ohmsweetohm1 · June 12, 2020, 3:46pm

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

What am I doing wrong?

klaff · June 12, 2020, 4:07pm

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

ohmsweetohm1 · June 12, 2020, 4:25pm

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.

klaff · June 12, 2020, 5:51pm

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.

klaff · June 12, 2020, 6:10pm

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]

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.

klaff · June 12, 2020, 6:20pm

ohmsweetohm1 · June 12, 2020, 6:28pm

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

klaff · June 12, 2020, 7:56pm

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

apo383 · June 12, 2020, 10:00pm

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.

ohmsweetohm1 · June 12, 2020, 10:58pm

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.

ChrisRackauckas · June 13, 2020, 8:24am

Hey,

When using the DEDataArray, you have to remember to update all of the cache variables. Notice the docs show the full_cache iterator.
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[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

ohmsweetohm1 · June 13, 2020, 9:24am

Thank you!

ohmsweetohm1 · June 14, 2020, 4:35pm

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

ChrisRackauckas · June 14, 2020, 4:47pm

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.