ContinuousCallback in DifferentialEquations.jl for the The Izhikevich Model

Specific Domains Modelling & Simulations
1

Greetings,
Is it possible to get the The Izhikevich Model to work with ContinuousCallback as the DiscreteCallback is not accurate near the zero-crossing.

# function reset!(integrator)
   # @show integrator.u[1], integrator.t
   # ...
(integrator.u[1], integrator.t) = (34.809473582040724, 2.596486286807452)
(integrator.u[1], integrator.t) = (32.229505727923915, 4.047660542741739)
(integrator.u[1], integrator.t) = (35.37791747113174, 5.655299412290314)
(integrator.u[1], integrator.t) = (36.029091214255416, 7.463878908313585)
(integrator.u[1], integrator.t) = (36.72604067495955, 9.567043050470012)
(integrator.u[1], integrator.t) = (31.099456800827046, 12.166997959461202)
...

I tried to change the parameters in this call:

cb =ContinuousCallback(thr, reset!;
    idxs = 1,
    interp_points = 2,
    abstol = eps(), reltol = 0, 
   )

but i always get:

┌ Warning: dt(5.684341886080802e-14) <= dtmin(5.684341886080802e-14) at t=2.8496719219500215, and step error estimate = 52.62708532815521. Aborting. There is either an error in your model specification or the true solution is unstable.

Here is the full code:

using DifferentialEquations
function izh!(du, u, p, t)
    a, b, c, d, I = p
    du[1] = 0.04 * u[1]^2 + 5.0 * u[1] + 140.0 - u[2] + I
    du[2] = a * (b * u[1] - u[2])

end
function thr(u, t, integrator)
    integrator.u[1] -30.0
end
function reset!(integrator)
    @show integrator.u[1], integrator.t
    integrator.u[1] = integrator.p[3]
    integrator.u[2] += integrator.p[4]
end
cb =ContinuousCallback(thr, reset!;
    idxs = 1,
    interp_points = 2,
    abstol = eps(), reltol = 0, 
   )
 solvers=[Tsit5()#= ,BS3() ,Rodas5(), ABDF2() ,  QNDF2() , Trapezoid() ,  TRBDF2()  , Rosenbrock23() =#]

for alg in solvers
    p = [0.02, 0.2, -50, 2.0, 10.0]
    u0 = [-60.0, p[2] * -60.0]
    tspan = (0.0, 300)
    prob = ODEProblem(izh!, u0, tspan, p , callback = cb)
    sol = solve(prob,alg,reltol=1e-3,abstol=1e-3);
end

many thanks!

2

I’m not sure exactly the answer, but invoking ModelingToolkit gives some magic that seems to work…

using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit: t_nounits as t, D_nounits as D
using Plots

pars = @parameters a=0.02 b=0.2 c=-50 d=2.0 I=10.0 
vars = @variables u1(t)=-60 u2(t)=-60*b

eqs = [
    D(u1) ~ 0.04 * u1^2 + 5.0 * u1 + 140.0 - u2 + I
    D(u2) ~ a * (b * u1 - u2)
]

root_eqs = [u1 ~ 30]
affect = [
    u1 ~ c
    u2 ~ Pre(u2) + d
]

@mtkcompile sys = System(eqs, t, vars, pars; continuous_events = root_eqs => affect)
prob = ODEProblem(sys, [], (0, 300))
sol = solve(prob) # Success
3

Some additional notes: I noticed you were using a loose tolerance of abstol, reltol = 1e-3, 1e-3. Using ModelingToolkitTolerances to run analysis shows this is likely to give high error/residual…

using ModelingToolkitTolerances
resids = analysis(prob, Tsit5())
plot(resids)

image
image597×403 18.1 KB

Running with tighter tolerances gives what seems to be a good result…

sol = solve(prob, Tsit5(); abstol=1e-9, reltol=1e-12)
plot(sol)

image
image618×397 41.8 KB

4

It’s a chaotic system so the analysis of ModelingToolkitTolerances won’t be so accurate as to stepper behavior. But…

You cannot resolve this to eps() after t=1 because machine epsilon is value dependent. So this is just not possible to achieve, and likely the reason why it’s going to not resolve.

5

Hi,
It turned out that i misued the integrator.u inside the thr function. If I use just u, it works fine. (not sure why)
Another question: the output plot does not contain data at the events

IzhikevichchModelTsit5_24nvautosave_
IzhikevichchModelTsit5_24nvautosave_600×400 50.4 KB

save_positions = (true, true) did not help. saveat works but I want to only save points at events.

Thanks again!

6

In the continuous callback? integrator.u is the value at integrator.t, where as u is the value at t. The latter is a constant so it’ll never have a root.

8

adding dense=false solved this issue (would be nicer if save_positions = (true, true) worked though)

IzhikevichchModelTsit5_25nvautosavel_
IzhikevichchModelTsit5_25nvautosavel_600×400 39.4 KB

Thank you all for your time

9

To save only at events, just set save_everystep=false and the default save_positions will be fine.