DiffEq tolerances on ContinuousCallbacks?

I’m troubleshooting some friction models and trying to understand what I’m seeing. I’ve made a MWE by taking the bouncing ball model from the DiffEq docs and changing the coefficient of restitution to 0.1 instead of 1.0 so that the bounces quickly get progressively smaller.

After the 5th bounce (with velocity down to 3e-4) the callback misses a zero cross and the ball falls through the earth. I know that if I want to handle a model that predicts an infinite series of smaller and smaller bounces, I will have to implement a heuristic to stop it before it fails due to numerical issues but I don’t know how to predict when this will fail, and it seems to be giving up earlier than necessary.

Is there a tolerance setting that controls the accuracy of the zero-cross finding process for the continuous callbacks?

MWE:

### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# ╔═╡ 2d729750-5455-11eb-323b-85261dc4f76a
begin # define independent package environment
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add("Plots"); using Plots
	Pkg.add("DifferentialEquations"); using DifferentialEquations
end

# ╔═╡ 59df1e80-5455-11eb-2323-1b07750a4d6b
function f(du,u,p,t)
  du[1] = u[2]
  du[2] = -p
end

# ╔═╡ 68dbca00-5455-11eb-36ba-ad0f2db6b034
function condition(u,t,integrator) # Event when event_f(u,t) == 0
  u[1]
end

# ╔═╡ 6ea96280-5455-11eb-07d4-adf2ce716897
function affect!(integrator)
  println("Bounce at t = ",integrator.t)
  integrator.u[2] = -0.1*integrator.u[2]
end

# ╔═╡ 7401cb00-5455-11eb-0b6d-973b7ed34eee
cb = ContinuousCallback(condition,affect!)

# ╔═╡ 7cdd3fc0-5455-11eb-15a2-07bf5923d0f7
function sim()
	u0 = [50.0,0.0]
	tspan = (0.0,5.0)
	p = 9.8
	prob = ODEProblem(f,u0,tspan,p)
	solve(prob,Tsit5(),callback=cb)
end

# ╔═╡ 2e6ff0ce-545a-11eb-3537-251e4eff7339
sol = sim()

# ╔═╡ 790adb00-5455-11eb-301c-61f7a5f0d3cf
function plotsol(tspan, title)
	p1 = plot(sol, vars=[1], tspan=tspan, 
		framestyle=:zerolines, ylabel="position", title=title)
	p2 = plot(sol, vars=[2], tspan=tspan,
		framestyle=:zerolines, ylabel="velocity")
	plot(p1,p2,layout=(2,1),link=:x, leg=false)
end

# ╔═╡ 63237f80-545b-11eb-37a5-830fefdc5c96
plotsol((0, 5), "can see first and second bounces clearly")

# ╔═╡ a830ce20-545b-11eb-3e3c-694bb6408666
plotsol((3.83,3.91), "last four bounce are in here")

# ╔═╡ bb97b18e-545b-11eb-0e1f-ebf59b70e16f
plotsol((3.895, 3.905), "last three bounces")

# ╔═╡ 596d1ea0-545c-11eb-07e3-3fcef590c7bd
plotsol((3.903, 3.905), "last two bounces")

# ╔═╡ cb3366c0-545c-11eb-2345-571ee3526ba5
plotsol((3.90415, 3.9043), "last bounce")

# ╔═╡ Cell order:
# ╠═59df1e80-5455-11eb-2323-1b07750a4d6b
# ╠═68dbca00-5455-11eb-36ba-ad0f2db6b034
# ╠═6ea96280-5455-11eb-07d4-adf2ce716897
# ╠═7401cb00-5455-11eb-0b6d-973b7ed34eee
# ╠═7cdd3fc0-5455-11eb-15a2-07bf5923d0f7
# ╠═63237f80-545b-11eb-37a5-830fefdc5c96
# ╠═a830ce20-545b-11eb-3e3c-694bb6408666
# ╠═bb97b18e-545b-11eb-0e1f-ebf59b70e16f
# ╠═596d1ea0-545c-11eb-07e3-3fcef590c7bd
# ╠═cb3366c0-545c-11eb-2345-571ee3526ba5
# ╠═2e6ff0ce-545a-11eb-3537-251e4eff7339
# ╠═790adb00-5455-11eb-301c-61f7a5f0d3cf
# ╠═2d729750-5455-11eb-323b-85261dc4f76a
2 Likes

Which version of DiffEqBase?

DiffEqBase v6.53.4

Julia Version 1.6.0-beta1.0
Commit b84990e1ac (2021-01-08 12:42 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-8650U CPU @ 1.90GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.0 (ORCJIT, skylake)
Environment:
  JULIA_NUM_THREADS = 8
  JULIA_REVISE_WORKER_ONLY = 1

Status `C:\Users\klaffedk\AppData\Local\Temp\jl_4bVRB6\Manifest.toml`
  [c3fe647b] AbstractAlgebra v0.12.0
  [79e6a3ab] Adapt v3.0.0
  [ec485272] ArnoldiMethod v0.0.4
  [4fba245c] ArrayInterface v2.14.13
  [4c555306] ArrayLayouts v0.4.12
  [aae01518] BandedMatrices v0.16.0
  [764a87c0] BoundaryValueDiffEq v2.7.0
  [fa961155] CEnum v0.4.1
  [d360d2e6] ChainRulesCore v0.9.24
  [35d6a980] ColorSchemes v3.10.2
  [3da002f7] ColorTypes v0.10.9
  [5ae59095] Colors v0.12.6
  [861a8166] Combinatorics v1.0.2
  [bbf7d656] CommonSubexpressions v0.3.0
  [34da2185] Compat v3.25.0
  [187b0558] ConstructionBase v1.0.0
  [d38c429a] Contour v0.5.7
  [adafc99b] CpuId v0.2.2
  [9a962f9c] DataAPI v1.4.0
  [864edb3b] DataStructures v0.18.8
  [e2d170a0] DataValueInterfaces v1.0.0
  [bcd4f6db] DelayDiffEq v5.28.1
  [2b5f629d] DiffEqBase v6.53.4
  [459566f4] DiffEqCallbacks v2.16.0
  [5a0ffddc] DiffEqFinancial v2.4.0
  [c894b116] DiffEqJump v6.12.2
  [77a26b50] DiffEqNoiseProcess v5.5.0
  [055956cb] DiffEqPhysics v3.8.0
  [163ba53b] DiffResults v1.0.3
  [b552c78f] DiffRules v1.0.2
  [0c46a032] DifferentialEquations v6.16.0
  [c619ae07] DimensionalPlotRecipes v1.2.0
  [b4f34e82] Distances v0.10.0
  [31c24e10] Distributions v0.24.10
  [ffbed154] DocStringExtensions v0.8.3
  [d4d017d3] ExponentialUtilities v1.8.0
  [e2ba6199] ExprTools v0.1.3
  [c87230d0] FFMPEG v0.4.0
  [9aa1b823] FastClosures v0.3.2
  [1a297f60] FillArrays v0.10.2
  [6a86dc24] FiniteDiff v2.7.2
  [53c48c17] FixedPointNumbers v0.8.4
  [59287772] Formatting v0.4.2
  [f6369f11] ForwardDiff v0.10.14
  [069b7b12] FunctionWrappers v1.1.1
  [28b8d3ca] GR v0.53.0
  [01680d73] GenericSVD v0.3.0
  [5c1252a2] GeometryBasics v0.3.5
  [42e2da0e] Grisu v1.0.0
  [cd3eb016] HTTP v0.8.19
  [615f187c] IfElse v0.1.0
  [d25df0c9] Inflate v0.1.2
  [83e8ac13] IniFile v0.5.0
  [c8e1da08] IterTools v1.3.0
  [42fd0dbc] IterativeSolvers v0.9.0
  [82899510] IteratorInterfaceExtensions v1.0.0
  [692b3bcd] JLLWrappers v1.2.0
  [682c06a0] JSON v0.21.1
  [b964fa9f] LaTeXStrings v1.2.0
  [2ee39098] LabelledArrays v1.4.0
  [23fbe1c1] Latexify v0.14.7
  [093fc24a] LightGraphs v1.3.4
  [d3d80556] LineSearches v7.1.1
  [bdcacae8] LoopVectorization v0.8.26
  [1914dd2f] MacroTools v0.5.6
  [739be429] MbedTLS v1.0.3
  [442fdcdd] Measures v0.3.1
  [e1d29d7a] Missings v0.4.4
  [961ee093] ModelingToolkit v4.5.0
  [46d2c3a1] MuladdMacro v0.2.2
  [f9640e96] MultiScaleArrays v1.8.1
  [d41bc354] NLSolversBase v7.7.1
  [2774e3e8] NLsolve v4.5.1
  [77ba4419] NaNMath v0.3.5
  [8913a72c] NonlinearSolve v0.3.4
  [6fe1bfb0] OffsetArrays v1.3.1
  [bac558e1] OrderedCollections v1.3.2
  [1dea7af3] OrdinaryDiffEq v5.49.1
  [90014a1f] PDMats v0.10.1
  [65888b18] ParameterizedFunctions v5.7.0
  [d96e819e] Parameters v0.12.1
  [69de0a69] Parsers v1.0.15
  [ccf2f8ad] PlotThemes v2.0.0
  [995b91a9] PlotUtils v1.0.10
  [91a5bcdd] Plots v1.10.1
  [7f904dfe] PlutoUI v0.6.11
  [e409e4f3] PoissonRandom v0.4.0
  [1fd47b50] QuadGK v2.4.1
  [74087812] Random123 v1.2.0
  [fb686558] RandomExtensions v0.4.3
  [e6cf234a] RandomNumbers v1.4.0
  [3cdcf5f2] RecipesBase v1.1.1
  [01d81517] RecipesPipeline v0.2.1
  [731186ca] RecursiveArrayTools v2.10.0
  [f2c3362d] RecursiveFactorization v0.1.6
  [189a3867] Reexport v0.2.0
  [ae029012] Requires v1.1.2
  [ae5879a3] ResettableStacks v1.1.0
  [79098fc4] Rmath v0.6.1
  [7e49a35a] RuntimeGeneratedFunctions v0.5.0
  [21efa798] SIMDPirates v0.8.26
  [476501e8] SLEEFPirates v0.5.5
  [1bc83da4] SafeTestsets v0.0.1
  [6c6a2e73] Scratch v1.0.3
  [efcf1570] Setfield v0.7.0
  [992d4aef] Showoff v0.3.2
  [699a6c99] SimpleTraits v0.9.3
  [a2af1166] SortingAlgorithms v0.3.1
  [47a9eef4] SparseDiffTools v1.12.0
  [276daf66] SpecialFunctions v1.2.1
  [90137ffa] StaticArrays v0.12.5
  [2913bbd2] StatsBase v0.33.2
  [4c63d2b9] StatsFuns v0.9.6
  [9672c7b4] SteadyStateDiffEq v1.6.0
  [789caeaf] StochasticDiffEq v6.30.1
  [09ab397b] StructArrays v0.4.2
  [c3572dad] Sundials v4.4.1
  [fd094767] Suppressor v0.2.0
  [d1185830] SymbolicUtils v0.6.3
  [3783bdb8] TableTraits v1.0.0
  [bd369af6] Tables v1.2.2
  [a759f4b9] TimerOutputs v0.5.7
  [a2a6695c] TreeViews v0.3.0
  [3a884ed6] UnPack v1.0.2
  [1986cc42] Unitful v1.5.0
  [3d5dd08c] VectorizationBase v0.12.33
  [19fa3120] VertexSafeGraphs v0.1.2
  [700de1a5] ZygoteRules v0.2.1
  [6e34b625] Bzip2_jll v1.0.6+5
  [83423d85] Cairo_jll v1.16.0+6
  [5ae413db] EarCut_jll v2.1.5+1
  [2e619515] Expat_jll v2.2.7+6
  [b22a6f82] FFMPEG_jll v4.3.1+4
  [a3f928ae] Fontconfig_jll v2.13.1+14
  [d7e528f0] FreeType2_jll v2.10.1+5
  [559328eb] FriBidi_jll v1.0.5+6
  [0656b61e] GLFW_jll v3.3.2+1
  [d2c73de3] GR_jll v0.53.0+0
  [78b55507] Gettext_jll v0.20.1+7
  [7746bdde] Glib_jll v2.59.0+4
  [aacddb02] JpegTurbo_jll v2.0.1+3
  [c1c5ebd0] LAME_jll v3.100.0+3
  [dd4b983a] LZO_jll v2.10.0+3
  [dd192d2f] LibVPX_jll v1.9.0+1
  [e9f186c6] Libffi_jll v3.2.1+4
  [d4300ac3] Libgcrypt_jll v1.8.5+4
  [7e76a0d4] Libglvnd_jll v1.3.0+3
  [7add5ba3] Libgpg_error_jll v1.36.0+3
  [94ce4f54] Libiconv_jll v1.16.0+7
  [4b2f31a3] Libmount_jll v2.34.0+3
  [89763e89] Libtiff_jll v4.1.0+2
  [38a345b3] Libuuid_jll v2.34.0+7
  [e7412a2a] Ogg_jll v1.3.4+2
  [458c3c95] OpenSSL_jll v1.1.1+6
  [efe28fd5] OpenSpecFun_jll v0.5.3+4
  [91d4177d] Opus_jll v1.3.1+3
  [2f80f16e] PCRE_jll v8.42.0+4
  [30392449] Pixman_jll v0.40.0+0
  [ede63266] Qt_jll v5.15.2+1
  [f50d1b31] Rmath_jll v0.2.2+1
  [fb77eaff] Sundials_jll v5.2.0+1
  [a2964d1f] Wayland_jll v1.17.0+4
  [2381bf8a] Wayland_protocols_jll v1.18.0+4
  [02c8fc9c] XML2_jll v2.9.10+3
  [aed1982a] XSLT_jll v1.1.33+4
  [4f6342f7] Xorg_libX11_jll v1.6.9+4
  [0c0b7dd1] Xorg_libXau_jll v1.0.9+4
  [935fb764] Xorg_libXcursor_jll v1.2.0+4
  [a3789734] Xorg_libXdmcp_jll v1.1.3+4
  [1082639a] Xorg_libXext_jll v1.3.4+4
  [d091e8ba] Xorg_libXfixes_jll v5.0.3+4
  [a51aa0fd] Xorg_libXi_jll v1.7.10+4
  [d1454406] Xorg_libXinerama_jll v1.1.4+4
  [ec84b674] Xorg_libXrandr_jll v1.5.2+4
  [ea2f1a96] Xorg_libXrender_jll v0.9.10+4
  [14d82f49] Xorg_libpthread_stubs_jll v0.1.0+3
  [c7cfdc94] Xorg_libxcb_jll v1.13.0+3
  [cc61e674] Xorg_libxkbfile_jll v1.1.0+4
  [12413925] Xorg_xcb_util_image_jll v0.4.0+1
  [2def613f] Xorg_xcb_util_jll v0.4.0+1
  [975044d2] Xorg_xcb_util_keysyms_jll v0.4.0+1
  [0d47668e] Xorg_xcb_util_renderutil_jll v0.3.9+1
  [c22f9ab0] Xorg_xcb_util_wm_jll v0.4.1+1
  [35661453] Xorg_xkbcomp_jll v1.4.2+4
  [33bec58e] Xorg_xkeyboard_config_jll v2.27.0+4
  [c5fb5394] Xorg_xtrans_jll v1.4.0+3
  [3161d3a3] Zstd_jll v1.4.5+2
  [0ac62f75] libass_jll v0.14.0+4
  [f638f0a6] libfdk_aac_jll v0.1.6+4
  [b53b4c65] libpng_jll v1.6.37+6
  [f27f6e37] libvorbis_jll v1.3.6+6
  [1270edf5] x264_jll v2020.7.14+2
  [dfaa095f] x265_jll v3.0.0+3
  [d8fb68d0] xkbcommon_jll v0.9.1+5
  [0dad84c5] ArgTools
  [56f22d72] Artifacts
  [2a0f44e3] Base64
  [ade2ca70] Dates
  [8bb1440f] DelimitedFiles
  [8ba89e20] Distributed
  [f43a241f] Downloads
  [9fa8497b] Future
  [b77e0a4c] InteractiveUtils
  [b27032c2] LibCURL
  [76f85450] LibGit2
  [8f399da3] Libdl
  [37e2e46d] LinearAlgebra
  [56ddb016] Logging
  [d6f4376e] Markdown
  [a63ad114] Mmap
  [ca575930] NetworkOptions
  [44cfe95a] Pkg
  [de0858da] Printf
  [3fa0cd96] REPL
  [9a3f8284] Random
  [ea8e919c] SHA
  [9e88b42a] Serialization
  [1a1011a3] SharedArrays
  [6462fe0b] Sockets
  [2f01184e] SparseArrays
  [10745b16] Statistics
  [4607b0f0] SuiteSparse
  [fa267f1f] TOML
  [a4e569a6] Tar
  [8dfed614] Test
  [cf7118a7] UUIDs
  [4ec0a83e] Unicode
  [e66e0078] CompilerSupportLibraries_jll
  [deac9b47] LibCURL_jll
  [29816b5a] LibSSH2_jll
  [c8ffd9c3] MbedTLS_jll
  [14a3606d] MozillaCACerts_jll
  [4536629a] OpenBLAS_jll
  [bea87d4a] SuiteSparse_jll
  [83775a58] Zlib_jll
  [8e850ede] nghttp2_jll

I found the abstol argument for the callback but even setting it ridiculously small (like 1e-40) doesn’t change the behavior. Also, I found this language interesting:

  • abstol=1e-14 & reltol=0: These are used to specify a tolerance from zero for the rootfinder:
    if the starting condition is less than the tolerance from zero, then no root will be detected.
    This is to stop repeat events happening just after a previously rootfound event. =#

Immediately after the callback, the position value is very close to zero (ideally it would be zero and the values around each of the bounces are reasonably close to eps(), but because the velocity is upwards, that position grows before falling back to zero. Here’s the plot of the last detected bounce and the next bounce which is not.

image

So if “starting condition” is the value of the condition expression immediately after callback execution, its value is irrelevant to how long it will be before the crossing will occur. Anyway I’ll look through the code in callbacks.jl a bit and see if anything jumps out at me.

Thanks.

Seems like a particularly hard case because of the very small time change associated with it. The typical issue here is the following. After an event, we want to make sure than manifold constraints are satisfied, i.e. if you want to inside of a zone you should stay inside. So we take the left side of a branching nonlinear solve, so it’s either exactly the zero or one floating point before. Given that, there’s no abstol anymore, and it should be right next to it.

But the hard part is then the next detection. If you didn’t get the zero exactly, are you positive or negative at the start point? Because you want to look for a sign change, you want a sign there. Normally we split the interval into 10 parts and the interval isn’t right at the start, but in this case, it’s right at the start. So you need a side on that start point. In the bouncing ball case, you want it to be positive. In something like a Poincure section, the user is just making a save at the zero and continuing, so there is no “bounce”, so right after the event, even if you pull to before the event, you want that to count as slightly negative (:scream:). So what we have to do is use the derivative at that starting point as our approximation for whether the condition is going up or down as our proxy for whether we should be treating it as slightly positive or slightly negative at the start.

It turns out that dt was too large in this case, using a point from where it has already gone negative, and boom. The fix is the following PR:

There’s a more interesting fix I was playing with, but I think it would need more playing to make it a reality:

Hopefully that’s enlightening as to why it’s such a complicated calculation though.

If you could contribute an event detection test DiffEqBase.jl/event_detection_tests.jl at master · SciML/DiffEqBase.jl · GitHub that would be great.

I was just reading through determine_event_occurance and coming to grips with the issues. Also that line with the semi-hardcoded “slightly in future” caught my eye and I was about to do some calculations on my solve points so I was pleased to learn I was on the right track when I saw your response.

Yes, I’ll add an event detection test.

This is easy enough to work around with a speed threshold which is practically a good idea anyway because real systems (at least the macroscopic ones I work with) don’t really keep bouncing but I’m struggling to think of a general way to pick that threshold. It’s essentially the same problem you have in trying to set a default absolute tolerance for general problems. It’ll be clearer in the morning.

Thanks again.

Real systems also aren’t such simple ODEs that dt keeps growing because the integrator can solve them without error (so adaptivity keeps going “oh hey, this is easy” while the events keep getting harder).

1 Like