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
). 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.