# DifferentialEquations error: "Event repeated at the same time"

Hi,

I’m trying to reproduce results presented in Roup et al. where a verge and foliot mechanism used in clocks is simulated. I tried to use a VectorContinuousCallback to determine the times at which collision events happen but I am getting an “Event repeated at the same time” error. This error occurs for a long enough integration time interval and/or a small enough timestep. Here’s the faulty code:

``````const torque = 1.0
const Ic = 10.0
const Iv = 0.15
const r_c = 1.0
const r_v = 0.3
const αc = 24.0 * pi/180
const e_r = 0.05
const nteeth = round(Int, 2pi/αc)
const Mv = Iv/r_v^2
const Mc = Ic/r_c^2
const Gc = Mv*(1 + e_r)/r_c/(Mv + Mc)
const Gv = Mc*(1 + e_r)/r_v/(Mv + Mc)
const αv = Gv*r_c/(1-e_r)*αc/2

p = (torque, Ic)

function f(du,u,p,t)
du[1] = u[3]
du[2] = u[4]
du[3] = p[1]/p[2] # p =(torque, crown inertia)
du[4] = 0.0
end

#out = vector of length 2nteeth
function condition(out,u,t,integrator)
for m in 1:nteeth
out[m] = r_c*sin(u[1] - (m-1)*αc) - r_v*tan(u[2] + αv/2)#upper
out[m + nteeth] = r_c*sin((m-1)*αc - u[1]) - r_v*tan(-u[2] + αv/2) #lower
end
end

function affect!(integrator, m)
c_up = r_c*integrator.u[3] - r_v*integrator.u[4] > 0.0 && ( (m-1)*αc <= mod(integrator.u[1] + αc/2, 2π) <= m*αc)#upper
c_dn = r_c*integrator.u[3] + r_v*integrator.u[4] > 0.0 && ( (m-1)*αc <= mod(integrator.u[1] + αc/2 + π, 2π) <= m*αc ) #lower
if m <= nteeth  && c_up
integrator.u[3] += -r_c*Gc*integrator.u[3] + r_v*Gc*integrator.u[4]
integrator.u[4] += + r_c*Gv*integrator.u[3] - r_v*Gv*integrator.u[4]
elseif m > nteeth && c_dn
integrator.u[3] +=  -r_c*Gc*integrator.u[3] - r_v*Gc*integrator.u[4]
integrator.u[4] +=  -r_c*Gv*integrator.u[3] - r_v*Gv*integrator.u[4]
end
end

cb = VectorContinuousCallback(condition,affect!,2nteeth)

u0 = [0, 0.0, 0.0, 0.0]
prob = ODEProblem(f,u0,(0.0,10.0),p)
sol = solve(prob, callback=cb, dt=1e-2, adaptive=false)
plot(sol,layout=(4,1))
``````

and the error:

``````Event repeated at the same time. Please report this error

Stacktrace:
[1] error(::String) at .\error.jl:33
[5] loopfooter! at C:\Users\bjd3\.julia\packages\OrdinaryDiffEq\kFkXd\src\integrators\integrator_utils.jl:166 [inlined]
[7] __solve(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Float64,Float64},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}}; kwargs::Base.Iterators.Pairs{Symbol,Any,NTuple{5,Symbol},NamedTuple{(:default_set, :second_time, :callback, :dt, :adaptive),Tuple{Bool,Bool,VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},Float64,Bool}}}) at C:\Users\bjd3\.julia\packages\OrdinaryDiffEq\kFkXd\src\solve.jl:5
[8] __solve(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Float64,Float64},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::Nothing; default_set::Bool, kwargs::Base.Iterators.Pairs{Symbol,Any,NTuple{4,Symbol},NamedTuple{(:second_time, :callback, :dt, :adaptive),Tuple{Bool,VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},Float64,Bool}}}) at C:\Users\bjd3\.julia\packages\DifferentialEquations\2DMrQ\src\default_solve.jl:7
[9] #__solve#467 at C:\Users\bjd3\.julia\packages\DiffEqBase\k3AhB\src\solve.jl:194 [inlined]
[10] (::DiffEqBase.var"#461#462"{ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Float64,Float64},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},Tuple{}})() at C:\Users\bjd3\.julia\packages\DiffEqBase\k3AhB\src\solve.jl:49
[11] maybe_with_logger at C:\Users\bjd3\.julia\packages\DiffEqBase\k3AhB\src\utils.jl:259 [inlined]
[12] solve_call(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Float64,Float64},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}; merge_callbacks::Bool, kwargs::Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol,Symbol},NamedTuple{(:callback, :dt, :adaptive),Tuple{VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},Float64,Bool}}}) at C:\Users\bjd3\.julia\packages\DiffEqBase\k3AhB\src\solve.jl:48
[13] #solve#463 at C:\Users\bjd3\.julia\packages\DiffEqBase\k3AhB\src\solve.jl:70 [inlined]
[14] top-level scope at In[8]:48
``````

I’m new to solving differential equations with events so I guess my first question is: did I properly set the problem up? In particular, the collision condition is determined by 3 sub-conditions, 2 booleans that I put in affect! and 1 that is root-findable that I put in the condition function. Is that the right way to deal with conditions that are of different nature?

If everything looks alright, could it be a bug in DifferentialEquations?
Thank you,

Bastien
P.S. Sorry this example isn’t very minimal, I didn’t want to remove parts of the models that are necessary to observe the expected periodic motion of the mechanism. Also, I didn’t write down the equations from the original article but I’d be happy to do so if anyone’s interested.

When I run this code, it seems to be working. It has been running for a few hours now, I’ll update if it finishes. The error message in you mentioned does not exist in current versions of `DifferentialEquations`, try updating your installed packages with `]update` in the REPL.

you can check the current versions with `]st` in the REPL:

``````(jl_KbfOVI) pkg> st
Status `/tmp/jl_KbfOVI/Project.toml`
[0c46a032] DifferentialEquations v6.15.0
[1dea7af3] OrdinaryDiffEq v5.46.0
``````

Since you specified that the callback is continuous (`VectorContinuousCallback`), `affect!` will be called for every zero-crossing in the output of `condition`. As long as every zero crossing is interpreted correctly by `affect!` and it should only be called on zero crossings, I think this callback structure should work fine.

There were some bugs like this, but they were fixed. Make sure you’re running an updated version.

Thanks @contradict and @ChrisRackauckas, I should have posted the versions. I was getting the error message in IJulia using DifferentialEquations v6.12.0.
The funny thing is that `Pkg.status()` points to a different path in IJulia and in the REPL, therefore pointing to different versions of `DifferentialEquations`. In the REPL, I got:

``````(@v1.5) pkg> st
Status `C:\Users\bjd3\.julia\environments\v1.5\Project.toml`
[0c46a032] DifferentialEquations v6.15.0
``````

whereas in IJulia:

``````Status `K:\bjd3\Project.toml`
[0c46a032] DifferentialEquations v6.12.0
``````

I updated to v6.15.0 in IJulia and it’s been running for some time now. I’ll update my post if the error message shows up again. @ChrisRackauckas is there a place in the doc that states which version fixes the bugs? I looked for a changelog but I couldn’t find it. It’d be great if there was a warning message in the section about callbacks that warn users of possible bugs prior to v6.15.0. That being said, `DifferentialEquations` is awesome!

This might not be the best place to do so but I wanted to give a big thanks to @ChrisRackauckas and all the contributors to `DifferentialEquations`: I am an experimental physicist with little background in numerics and switching to julia from matlab and python wasn’t an easy call at first. It was totally worth it , I got speedups in the 50-100X range that simply enabled work to be done. Maybe even more importantly, Chris’ videos are crystal clear and brought concepts like non-allocation to my attention. I learnt more in a month using julia than in a year using matlab…

2 Likes

The REPL is using the default Julia environment. You can read more about the idea of an environment here: https://docs.julialang.org/en/v1/manual/code-loading/#Environments

When you start a new IJulia notebook, it runs Julia in an environment with the same directory as the notebook: https://julialang.github.io/IJulia.jl/stable/manual/usage/#Julia-projects

In any case, it sounds like you figured out how to get it updated!

This is a fun problem, I enjoyed looking over that paper.

Thanks for the links, I didn’t pay attention to that “detail” when I installed IJulia. My bad…

Did you let the simulations run by any chance? I’ve been running it for a couple hours on my computer and it’s still working on it. I will benchmark it a little more carefully when I have a second but it’s weird that it’s taking that long because the integration was done in a matter of a few minutes or less using v6.12.0 (when using a timespan and a timestep that worked).

1 Like

Just checked on it, still running. I’m a little surprised too, I figured it would have finished overnight. With only 1e3 timesteps to compute, something must be funny with the callback but I don’t see what.

What’s your final code here?

my code is currently:

``````const torque = 1.0
const Ic = 10.0
const Iv = 0.15
const r_c = 1.0
const r_v = 0.3
const αc = 24.0 * pi/180
const e_r = 0.05
const nteeth = round(Int, 2pi/αc)
const Mv = Iv/r_v^2
const Mc = Ic/r_c^2
const Gc = Mv*(1 + e_r)/r_c/(Mv + Mc)
const Gv = Mc*(1 + e_r)/r_v/(Mv + Mc)
const αv = Gv*r_c/(1-e_r)*αc/2

p = (torque, Ic)

function f(du,u,p,t)
du[1] = u[3]
du[2] = u[4]
du[3] = p[1]/p[2] # p =(torque, crown inertia)
du[4] = 0.0
end

#out = vector of length 2nteeth
function condition(out,u,t,integrator)
for m in 1:nteeth
out[m] = r_c*sin(u[1] - (m-1)*αc) - r_v*tan(u[2] + αv/2)#upper
out[m + nteeth] = r_c*sin((m-1)*αc - u[1]) - r_v*tan(-u[2] + αv/2) #lower
end
end

function affect!(integrator, m)
c_up = r_c*integrator.u[3] - r_v*integrator.u[4] > 0.0 && ( (m-1)*αc <= mod(integrator.u[1] + αc/2, 2π) <= m*αc)#upper
c_dn = r_c*integrator.u[3] + r_v*integrator.u[4] > 0.0 && ( (m-1)*αc <= mod(integrator.u[1] + αc/2 + π, 2π) <= m*αc ) #lower
if  m <= nteeth  && c_up
integrator.u[3] += -r_c*Gc*integrator.u[3] + r_v*Gc*integrator.u[4]
integrator.u[4] += + r_c*Gv*integrator.u[3] - r_v*Gv*integrator.u[4]
elseif m > nteeth && c_dn
integrator.u[3] +=  -r_c*Gc*integrator.u[3] - r_v*Gc*integrator.u[4]
integrator.u[4] +=  -r_c*Gv*integrator.u[3] - r_v*Gv*integrator.u[4]
end
end

cb = VectorContinuousCallback(condition,affect!,2nteeth)

u0 = [0, 0.0, 0.0, 3.0]
prob = ODEProblem(f,u0,(0.0,2.0),p)
sol = solve(prob, callback=cb) #, adaptive=false
plot(sol,layout=(4,1))
``````

The problem is that it takes several hours to run (I didn’t time it properly, will do today) and exit with a maxiter error. Compared to the previously posted code, I reduced the timespan and made the solver adaptative. Concerning the adaptiveness, I hoped that the solver would be smart enough to go fast when there’s no discontinuity but I ended up with a sol with 1999998-elements (hence the maxiter that stops the integration at t=0.4 instead of 2). In the doc examples, the bouncing ball has an adaptative solver but the bouncing ball with walls hasn’t but there isn’t a clear explanation why. My plan was to look into the faq when I have some time to try to speed it up.

I’m bookmarking this to look at it after grades are in. It might be getting trapped somewhere? Do `@show t` in your `f` function and see how time progresses.

It definitely looks like the solver gets trapped: when I used the code posted above, I got this solution that definitely looks stuck at t=0.413367…

``````retcode: MaxIters
Interpolation: automatic order switching interpolation
t: 1999998-element Array{Float64,1}:
0.0
0.00033303728989641926
0.003663410188860612
0.036967139178502535
0.09147424476687356
0.09147424476687356
0.11091476940398995
0.11091476940398995
0.20226200532656313
0.20226200532656313
0.2554499762884815
0.2554499762884815
0.2858965285936999
⋮
0.413367402556334
0.413367402556334
0.413367402556334
0.413367402556334
0.413367402556334
0.413367402556334
0.413367402556334
0.413367402556334
0.413367402556334
0.413367402556334
0.413367402556334
0.413367402556334
``````

and using `@show t` in `f`, we see that the integration progresses well at first but eventually gets stuck around 0.4133…

``````t = 0.0
t = 0.0161
t = 0.0327
t = 0.09000000000000001
t = 0.09800255409045097
t = 0.1
t = 0.1
t = 0.014727353407466691
t = 0.02991207803876775
t = 0.08232682029018647
t = 0.08964709620648706
t = 0.09147424476687385
t = 0.09147424476687385
t = 0.014727353407466691
t = 0.02991207803876775
t = 0.08232682029018647
t = 0.08964709620648706
t = 0.09147424476687385
t = 0.09147424476687385
t = 0.09147424476687385
t = 0.10757424476687386
t = 0.12417424476687386
t = 0.18147424476687385
t = 0.18947679885732482
t = 0.19147424476687386
t = 0.19147424476687386
t = 0.0946041692334496
t = 0.09783129632321091
t = 0.10897071694027863
t = 0.1105264554398313
t = 0.11091476940399027
t = 0.11091476940399027
t = 0.0946041692334496
t = 0.09783129632321091
t = 0.10897071694027863
t = 0.1105264554398313
t = 0.11091476940399027
t = 0.11091476940399027
t = 0.11091476940399027
t = 0.12701476940399026
t = 0.14361476940399026
t = 0.20091476940399028
t = 0.20891732349444125
t = 0.2109147694039903
t = 0.2109147694039903
t = 0.12562167438752458
t = 0.14078531555067175
t = 0.1931272817343063
t = 0.2004373936991421
t = 0.20226200532656363
t = 0.20226200532656363
t = 0.12562167438752458
t = 0.14078531555067175
t = 0.1931272817343063
t = 0.2004373936991421
t = 0.20226200532656363
t = 0.20226200532656363
t = 0.20226200532656363
t = 0.21836200532656364
t = 0.23496200532656364
t = 0.29226200532656366
t = 0.3002645594170146
t = 0.30226200532656367
t = 0.30226200532656367
t = 0.21082526865143253
t = 0.21965447183111103
t = 0.2501311791922904
t = 0.2543875753381313
t = 0.25544997628848226
t = 0.25544997628848226
t = 0.21082526865143253
t = 0.21965447183111103
t = 0.2501311791922904
t = 0.2543875753381313
t = 0.25544997628848226
t = 0.25544997628848226
t = 0.25544997628848226
t = 0.27154997628848226
t = 0.28814997628848227
t = 0.3454499762884823
t = 0.3534525303789332
t = 0.3554499762884823
t = 0.3554499762884823
t = 0.26035187120962244
t = 0.26540599889228866
t = 0.28285187336317885
t = 0.2852883751800814
t = 0.2858965285937007
t = 0.2858965285937007
t = 0.26035187120962244
t = 0.26540599889228866
t = 0.28285187336317885
t = 0.2852883751800814
t = 0.2858965285937007
t = 0.2858965285937007
t = 0.2858965285937007
t = 0.3019965285937007
t = 0.3185965285937007
t = 0.3758965285937007
t = 0.3838990826841516
t = 0.3858965285937007
t = 0.3858965285937007
t = 0.2884789944932398
t = 0.29114166119711243
t = 0.30033267337373304
t = 0.30161629592328626
t = 0.3019366894604033
t = 0.3019366894604033
t = 0.2884789944932398
t = 0.29114166119711243
t = 0.30033267337373304
t = 0.30161629592328626
t = 0.3019366894604033
t = 0.3019366894604033
t = 0.3019366894604033
t = 0.3180366894604033
t = 0.3346366894604033
t = 0.39193668946040333
t = 0.3999392435508543
t = 0.4019366894604033
t = 0.4019366894604033
t = 0.30302852607470415
t = 0.30415427065839323
t = 0.3080401239502838
t = 0.30858282466853704
t = 0.3087182833380483
t = 0.3087182833380483
t = 0.30302852607470415
t = 0.30415427065839323
t = 0.3080401239502838
t = 0.30858282466853704
t = 0.3087182833380483
t = 0.3087182833380483
t = 0.3087182833380483
t = 0.3248182833380483
t = 0.3414182833380483
t = 0.3987182833380483
t = 0.4067208374284993
t = 0.40871828333804827
t = 0.40871828333804827
t = 0.30919250202321813
t = 0.3096814479967473
t = 0.31136919524272444
t = 0.31160490708612654
t = 0.3116637410099107
t = 0.3116637410099107
t = 0.30919250202321813
t = 0.3096814479967473
t = 0.31136919524272444
t = 0.31160490708612654
t = 0.3116637410099107
t = 0.3116637410099107
t = 0.3116637410099107
t = 0.3277637410099107
t = 0.3443637410099107
t = 0.4016637410099107
t = 0.4096662951003617
t = 0.41166374100991066
t = 0.41166374100991066
t = 0.31230724243837865
t = 0.3129707283832587
t = 0.31526095396407955
t = 0.31558080831119745
t = 0.31566064429232055
t = 0.31566064429232055
t = 0.31230724243837865
t = 0.3129707283832587
t = 0.31526095396407955
t = 0.31558080831119745
t = 0.31566064429232055
t = 0.31566064429232055
t = 0.31566064429232055
t = 0.33176064429232055
t = 0.34836064429232055
t = 0.40566064429232057
t = 0.41366319838277155
t = 0.4156606442923205
t = 0.4156606442923205
t = 0.33139144072283033
t = 0.3476107712039771
t = 0.4035967734069716
t = 0.41141581373747016
t = 0.41336745441971057
t = 0.41336745441971057
t = 0.33139144072283033
t = 0.3476107712039771
t = 0.4035967734069716
t = 0.41141581373747016
t = 0.41336745441971057
t = 0.41336745441971057
t = 0.41336745441971057
``````

where I used `sol = solve(prob, callback=cb, abstol=1e-3, reltol=1e-3, dt = 1e-1, adaptive=false)` . Note that I truncated the output of `@show` because you can already see the pattern forming but eventually t gets trapped on t = 0.41336745441971057.

Open an issue and I’ll investigate more. @kanav99 might want to take a look too. It might be repeatedly catching the same event, but not throwing the error? That would be odd.

I just posted the issue on github: https://github.com/SciML/DifferentialEquations.jl/issues/703

Thanks for looking at it.