DifferentialEquations error: "Event repeated at the same time"

Specific Domains Modelling & Simulations
1

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
 [2] apply_callback!(::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}},true,Array{Float64,1},Nothing,Float64,Tuple{Float64,Float64},Float64,Float64,Float64,Array{Array{Float64,1},1},OrdinaryDiffEq.ODECompositeSolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},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}},OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23Tableau{Float64},DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Tuple{Float64,Float64}},DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Tuple{Float64,Float64}},DefaultLinSolve,FiniteDiff.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},UnitRange{Int64},Nothing,Val{:forward},Float64},FiniteDiff.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:forward},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}}},DiffEqBase.DEStats},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23Tableau{Float64},DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Tuple{Float64,Float64}},DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Tuple{Float64,Float64}},DefaultLinSolve,FiniteDiff.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},UnitRange{Int64},Nothing,Val{:forward},Float64},FiniteDiff.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:forward},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(LinearAlgebra.opnorm),CallbackSet{Tuple{VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Float64,DiffEqBase.CallbackCache{Array{Float64,1},Array{Float64,1}},OrdinaryDiffEq.DefaultInit}, ::VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}, ::Float64, ::Float64, ::Int64) at C:\Users\bjd3\.julia\packages\DiffEqBase\k3AhB\src\callbacks.jl:758
 [3] handle_callbacks!(::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}},true,Array{Float64,1},Nothing,Float64,Tuple{Float64,Float64},Float64,Float64,Float64,Array{Array{Float64,1},1},OrdinaryDiffEq.ODECompositeSolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},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}},OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23Tableau{Float64},DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Tuple{Float64,Float64}},DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Tuple{Float64,Float64}},DefaultLinSolve,FiniteDiff.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},UnitRange{Int64},Nothing,Val{:forward},Float64},FiniteDiff.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:forward},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}}},DiffEqBase.DEStats},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23Tableau{Float64},DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Tuple{Float64,Float64}},DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Tuple{Float64,Float64}},DefaultLinSolve,FiniteDiff.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},UnitRange{Int64},Nothing,Val{:forward},Float64},FiniteDiff.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:forward},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(LinearAlgebra.opnorm),CallbackSet{Tuple{VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Float64,DiffEqBase.CallbackCache{Array{Float64,1},Array{Float64,1}},OrdinaryDiffEq.DefaultInit}) at C:\Users\bjd3\.julia\packages\OrdinaryDiffEq\kFkXd\src\integrators\integrator_utils.jl:252
 [4] _loopfooter!(::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}},true,Array{Float64,1},Nothing,Float64,Tuple{Float64,Float64},Float64,Float64,Float64,Array{Array{Float64,1},1},OrdinaryDiffEq.ODECompositeSolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},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}},OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23Tableau{Float64},DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Tuple{Float64,Float64}},DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Tuple{Float64,Float64}},DefaultLinSolve,FiniteDiff.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},UnitRange{Int64},Nothing,Val{:forward},Float64},FiniteDiff.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:forward},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}}},DiffEqBase.DEStats},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23Tableau{Float64},DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Tuple{Float64,Float64}},DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Tuple{Float64,Float64}},DefaultLinSolve,FiniteDiff.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},UnitRange{Int64},Nothing,Val{:forward},Float64},FiniteDiff.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:forward},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(LinearAlgebra.opnorm),CallbackSet{Tuple{VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Float64,DiffEqBase.CallbackCache{Array{Float64,1},Array{Float64,1}},OrdinaryDiffEq.DefaultInit}) at C:\Users\bjd3\.julia\packages\OrdinaryDiffEq\kFkXd\src\integrators\integrator_utils.jl:220
 [5] loopfooter! at C:\Users\bjd3\.julia\packages\OrdinaryDiffEq\kFkXd\src\integrators\integrator_utils.jl:166 [inlined]
 [6] solve!(::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}},true,Array{Float64,1},Nothing,Float64,Tuple{Float64,Float64},Float64,Float64,Float64,Array{Array{Float64,1},1},OrdinaryDiffEq.ODECompositeSolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},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}},OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23Tableau{Float64},DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Tuple{Float64,Float64}},DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Tuple{Float64,Float64}},DefaultLinSolve,FiniteDiff.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},UnitRange{Int64},Nothing,Val{:forward},Float64},FiniteDiff.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:forward},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}}},DiffEqBase.DEStats},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23Tableau{Float64},DiffEqBase.TimeGradientWrapper{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},Tuple{Float64,Float64}},DiffEqBase.UJacobianWrapper{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Tuple{Float64,Float64}},DefaultLinSolve,FiniteDiff.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},UnitRange{Int64},Nothing,Val{:forward},Float64},FiniteDiff.GradientCache{Nothing,Array{Float64,1},Array{Float64,1},Val{:forward},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(LinearAlgebra.opnorm),CallbackSet{Tuple{VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Float64,DiffEqBase.CallbackCache{Array{Float64,1},Array{Float64,1}},OrdinaryDiffEq.DefaultInit}) at C:\Users\bjd3\.julia\packages\OrdinaryDiffEq\kFkXd\src\solve.jl:425
 [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.

2

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.

3

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

4

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
5

The REPL is using the default Julia environment. You can read more about the idea of an environment here: Code Loading · The Julia Language

When you start a new IJulia notebook, it runs Julia in an environment with the same directory as the notebook: Using IJulia · IJulia

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.

6

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
7

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.

8

What’s your final code here?

9

Hi @ChrisRackauckas,

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.

10

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.

11

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.

12

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.

13

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

Thanks for looking at it.