Stepping on exactly the event point

[anirudhpammi]

I tried the code given on the link (https://docs.juliadiffeq.org/latest/features/callback_functions/). The Tsit5() adaptive algorithm doesn't work for this example. Can you please tell me why that is the case? When I disable the adaptive step size, specify the dt=1e-3 and solve this example, I get matching results.

Unrelated question: I have a different case, I need to check if one of the state variables of my system have reached a certain value. If I use the DiscreteCallback, I notice that there is an overshoot in my simulations. Should I be using the ContinuousCallback for better results?

[ChrisRackauckas]

I tried the code given on the link (https://docs.juliadiffeq.org/latest/features/callback_functions/). The Tsit5() adaptive algorithm doesn't work for this example. Can you please tell me why that is the case? When I disable the adaptive step size, specify the dt=1e-3 and solve this example, I get matching results.

Which example? Maybe @kanav99 could take a look.

If I use the DiscreteCallback, I notice that there is an overshoot in my simulations. Should I be using the ContinuousCallback for better results?

They do very different things. DiscreteCallbacks act at the end of a step, ContinuousCallback acts at a rootfinding condition. If you need to hit a specific value, like the bouncing ball, use the latter. If you want to do something like plot after each step, use the former.

[kanav99]

Hi! I just tried the example on OrdinaryDiffEq@master and DiffEqBase@master and it works for me. Did you make any changes to the code to cause error? Also are you on updated DifferentialEquations.jl?

[anirudhpammi]

@ChrisRackauckas : Thank you! That makes it clear. @kanav99 : It is the example of the bouncing ball, I am using the v6.9.0 of DifferentialEquations. The only change I made was to use a different plotting backend. I use the PyPlot backend and pygui option to plot the figures in an external window. If I use the Plots package and plot the solution using the command "plot(sol)" it seems to work fine. I guess it means that the plot command makes sense of the solution object and plots it in a meaningful way. Do you know if there is a way around this so that I can use just the sol.t and sol.u points to plot a meaningful answer without any backend magic? Thank you.

[anirudhpammi]

@ChrisRackauckas : Just one more follow up question, for the continuouscallback, what sort of root finding method do you use? Thanks a lot!

[kanav99]

That seems to be some plotting issue then. Let me check it. We use Bisection as of now for rootfinding.

[ChrisRackauckas]

You have to make sure you plot the discontinuity points on both sides to resolve it correctly. What does your plot look like? Sounds like you just plotted a standard plot over time, so you got a straight line between each point.

[anirudhpammi]

@ChrisRackauckas : Yes, that is exactly what I was doing and I get a straight line. I will look into it . I will read the documentation about the discontinuity points.

@kanav99 : Thank you, I seem to have a problem with the algorithm, at some time step my state variable steps exactly on the point which is a zero of the function it is supposed to minimize, and the bisection method fails to identity it as a solution. The documentation says that if the starting condition is < 10*eps() (by default) the root can't be found and since my starting point is exactly the root and event is not triggered. Is this intended? Could you suggest a workaround? Thanks!

[ChrisRackauckas]

Sounds like it was just a user plotting issue. Please feel free to discuss more if there's more to it!

[anirudhpammi]

@ChrisRackauckas : Sorry for the delay in replying. I did have a second question about the root finding algorithm: at some time step my state variable steps exactly on the point which is a zero of the function it is supposed to minimize, and the bisection method fails to identity it as a solution. The documentation says that if the starting condition is < 10*eps() (by default) the root can't be found and since my starting point is exactly the root and event is not triggered. Is this intended? Could you suggest a workaround? Thanks!

[ChrisRackauckas]

The start is treated differently, since otherwise a lot of callbacks would have weird side effects. For that reason we put in an initialize! option in the callback interface. If you want the effect to apply at the starting point, then you can add initialize = integrator -> affect!(integrator). By default it's a function that does nothing.

[anirudhpammi]

My starting point isn't the problem. The state variable exactly steps on the value that minimizes the function during the simulation. Since the step is exactly on the value, the root finding algorithm fails. I encountered a similar problem with a python package as well. Is that normal?

[ChrisRackauckas]

Hmm, that is odd. Can you share an example of that? We can probably find a way to fix it.

[anirudhpammi]

The brief description of the problem is in the pdf: Solver_stepping_exactly_on_root (1).pdf

using PyPlot, DifferentialEquations, BenchmarkTools
pygui(true)
close("all")

# Parameters
g = -0.4
α = 1
t_exc_1 = []
t_exc_2 = []

function Es(t,t_exc)
    t_redef = (-t_exc .+ t) .* Int.(t .>= t_exc)
    return ( length(t_redef) == 0 ? 0 : sum((g*α^2) .* exp.(-α .* t_redef) .* t_redef)  )
end

function neuron_activation!(du,u,p,t)
    # du -> rate of change, u -> state variable/initial conditions, p -> parameters, t -> time
    X = p[1]
    du[1] = X - u[1] + Es(t,t_exc_2)
    du[2] = X - u[2] + Es(t,t_exc_1)
end

function condition(u,t,integrator)
    (u[1]-1)
end

function affect!(integrator)
    integrator.u[1] = 0
    push!(t_exc_1,integrator.t)
end

function condition2(u,t,integrator)
    (u[2] - 1)
end

function affect2!(integrator)
    integrator.u[2] = 0
    push!(t_exc_2,integrator.t)
end

cb1 = ContinuousCallback(condition,affect!)
cb2 = ContinuousCallback(condition2,affect2!)

u0 = [0.0;0.001]
tspan = (0.0,200)
X_thresh = 1.3
p = [X_thresh]
prob = ODEProblem(neuron_activation!,u0,tspan,p);
alg = Vern7();
sol = solve(prob,alg,callback=CallbackSet(cb2,cb1),reltol=1e-6,abstol=1e-7)

figure()
plot(sol.t,sol.u)
T_p = diff(t_exc_1)[end]

[anirudhpammi]

@ChrisRackauckas : Did you have a chance to look at it? Thanks!

[ChrisRackauckas]

I don't understand the issue. I ran:

using DifferentialEquations, BenchmarkTools

# Parameters
g = -0.4
α = 1
t_exc_1 = []
t_exc_2 = []

function Es(t,t_exc)
    t_redef = (-t_exc .+ t) .* Int.(t .>= t_exc)
    return ( length(t_redef) == 0 ? 0 : sum((g*α^2) .* exp.(-α .* t_redef) .* t_redef)  )
end

function neuron_activation!(du,u,p,t)
    # du -> rate of change, u -> state variable/initial conditions, p -> parameters, t -> time
    X = p[1]
    du[1] = X - u[1] + Es(t,t_exc_2)
    du[2] = X - u[2] + Es(t,t_exc_1)
end

function condition(u,t,integrator)
    (u[1]-1)
end

function affect!(integrator)
    integrator.u[1] = 0
    @show integrator.t
    push!(t_exc_1,integrator.t)
end

function condition2(u,t,integrator)
    (u[2] - 1)
end

function affect2!(integrator)
    integrator.u[2] = 0
    @show integrator.t
    push!(t_exc_2,integrator.t)
end

cb1 = ContinuousCallback(condition,affect!)
cb2 = ContinuousCallback(condition2,affect2!)

u0 = [0.0;0.001]
tspan = (0.0,200)
X_thresh = 1.3
p = [X_thresh]
prob = ODEProblem(neuron_activation!,u0,tspan,p);
alg = Vern7();
sol = solve(prob,alg,callback=CallbackSet(cb2,cb1),reltol=1e-6,abstol=1e-7)

using Plots
Plots.plot(sol,plotdensity=10000)
T_p = diff(t_exc_1)[end]

and it seems both are hit around every 4 seconds?

1.4655846249344873
integrator.t = 1.4663557393427966
integrator.t = 4.833495929990829
integrator.t = 4.8337186499965545
integrator.t = 8.211811740475266
integrator.t = 8.21240923647408
integrator.t = 11.732253853253649
integrator.t = 11.732380259292999
integrator.t = 15.364787073434572
integrator.t = 15.364795807589566
integrator.t = 18.644838417286174
integrator.t = 18.644858635980874
integrator.t = 22.291288930455142
integrator.t = 22.291295362055575
integrator.t = 25.685231382573736
integrator.t = 25.685235779396574
integrator.t = 29.17958109620059
integrator.t = 29.17958169387737
integrator.t = 32.823674570321764
integrator.t = 32.82367534448865
integrator.t = 36.089201635021496
integrator.t = 36.08920181262052
integrator.t = 39.728371029287764
integrator.t = 39.72837107358062
integrator.t = 43.16127715469289
integrator.t = 43.16127720291633
integrator.t = 46.62742473319724
integrator.t = 46.62742473750797
integrator.t = 50.27959576551974
integrator.t = 50.2795957671167
integrator.t = 53.538355348904325
integrator.t = 53.538355349666766
integrator.t = 57.16759778651662
integrator.t = 57.167597786706295
integrator.t = 60.64257505981047
integrator.t = 60.642575060043335
integrator.t = 64.0778526434133
integrator.t = 64.07785264341841
integrator.t = 67.73482900472716
integrator.t = 67.73482900475796
integrator.t = 70.99237646760817
integrator.t = 70.99237646762145
integrator.t = 74.60900481375637
integrator.t = 74.60900481375741
integrator.t = 78.12344130481995
integrator.t = 78.12344130482177
integrator.t = 81.52765431951762
integrator.t = 81.52765431951794
integrator.t = 85.18679436444923
integrator.t = 85.18679436444924
integrator.t = 88.45232341879465
integrator.t = 88.45232341879466
integrator.t = 92.05297251007462
integrator.t = 92.05297251007464
integrator.t = 95.60380586730281
integrator.t = 95.60380586730284
integrator.t = 98.97691545616867
integrator.t = 98.97691545616868
integrator.t = 102.6358584814344
integrator.t = 102.63585848143441
integrator.t = 105.91899192836503
integrator.t = 105.91899192836505
integrator.t = 109.49992210691605
integrator.t = 109.49992210691606
integrator.t = 113.08169124894346
integrator.t = 113.08169124894347
integrator.t = 116.4246416280805
integrator.t = 116.4246416280805
integrator.t = 120.08152614993232
integrator.t = 120.08152614993232
integrator.t = 123.38982987810341
integrator.t = 123.38982987810341
integrator.t = 126.94812654085968
integrator.t = 126.94812654085968
integrator.t = 130.5546587497491
integrator.t = 130.5546587497491
integrator.t = 133.86940629010127
integrator.t = 133.86940629010127
integrator.t = 137.5229006709948
integrator.t = 137.5229006709948
integrator.t = 140.86345080246676
integrator.t = 140.86345080246676
integrator.t = 144.3968704210837
integrator.t = 144.3968704210837
integrator.t = 148.02234456114817
integrator.t = 148.02234456114817
integrator.t = 151.3118128187459
integrator.t = 151.3118128187459
integrator.t = 154.96087467905457
integrator.t = 154.96087467905457
integrator.t = 158.3369251574098
integrator.t = 158.3369251574098
integrator.t = 161.84410625747992
integrator.t = 161.84410625747992
integrator.t = 165.48322813590906
integrator.t = 165.48322813590906
integrator.t = 168.755097626517
integrator.t = 168.755097626517
integrator.t = 172.39787010241372
integrator.t = 172.39787010241372
integrator.t = 175.81174903239034
integrator.t = 175.81174903239034
integrator.t = 179.2915042355064
integrator.t = 179.29150423550644
integrator.t = 182.9401165481654
integrator.t = 182.9401165481654
integrator.t = 186.20046350416428
integrator.t = 186.20046350416428
integrator.t = 189.8350449767773
integrator.t = 189.8350449767773
integrator.t = 193.28922759355197
integrator.t = 193.28922759355197
integrator.t = 196.73990263202768
integrator.t = 196.73990263202768

[anirudhpammi]

Could you please try with u0 = [0.0;0.0] instead of u0=[0.0,0.001]. I notice that with the default value of reltol and abstol, it seems to work fine, but with the values I gave (reltol=1e-6,abstol=1e-7), the solutions doesn't work correctly.

[kanav99]

@anirudhpammi does this plot look good? I think this issue is fixed due to recent changes.

[ChrisRackauckas]

That looks correct to me.

[anirudhpammi]

Yes, it looks correct. thanks!