# DifferentialEquations; how to simulate limit stops

I’m trying to simulate a spring mass system in which the mass is constrained by limits stops or walls.

There’s a driving force which, if it pushes the mass to one wall or the other, essentially halts the system until the net force will take the mass off the wall. For MWE purposes, imagine a sinusoidal source of force which is of large enough magnitude to alternately pin the mass against one wall or the other, and the interesting behavior is what happens in between.

I’m having trouble figuring out how to arrange callbacks to achieve this, partly because the condition to get off a wall is not dependent on the state (which in the MWE would only include position and velocity, but not the force acting on the mass).

Any thoughts about how to structure such a problem?

``````using DifferentialEquations, Plots
plotly()

fdrive(t) = t%4.0>2.0 ? 5.0 : -5.0

function spring_mass_ODE!(du,u,p,t)
ẋ, x = u
mass, k, kd = p
du[1] = 1/mass*(-k*x-kd*ẋ+fdrive(t))
du[2] = ẋ
end

function smsim()
tspan = [0.0, 10.0]
u0 = [1.0, 1.0]
p=[1.0, 1.0, 0.2]
prob = ODEProblem(spring_mass_ODE!, u0, tspan,p)
sol = solve(prob,Tsit5(), abstol=1e-7,reltol=1e-7)
end

function plotsmsim(sol)
p1=plot(sol, vars=[1,],ylabel="velocity")
p2=plot(sol, vars=[2,],ylabel="position")
end
``````

If you run `smsim() |> plotsmsim` you see a position that ranges from -6.something to over 6.something. I want to place walls at say, +/- 3, so that it runs against the wall at -3, then once the force reverses swings up to +3, sticks there a while, and so on.

What do you mean? So have a callback that starts the stick, and have a callback that unsticks using? A timer? You just need to specify what that model is.

One easy way that doesn’t need callback is to put a very stiff spring at each of your limits:

``````# useful for defining piecewise functions, 0 for negative argument, 1 for positive
heaviside(x) = 0.5 * (sign(x) + 1)

function spring_mass_ODE!(du,u,p,t)
ẋ, x = u
mass, k, kd, stop_rate = p
negative_stop_force =  heaviside(-3 - x) * stop_rate
positive_stop_force = -heaviside(x - 3) * stop_rate
du[1] = 1/mass*(-k*x-kd*ẋ+negative_stop_force+positive_stop_force+fdrive(t))
du[2] = ẋ
end

function smsim()
tspan = [0.0, 10.0]
u0 = [1.0, 1.0]
p=[1.0, 1.0, 0.2, 100.0] # new fourth parameter for stop spring strength
prob = ODEProblem(spring_mass_ODE!, u0, tspan,p)
sol = solve(prob,Tsit5(), abstol=1e-7,reltol=1e-7)
end
``````

Another solution that would need callbacks is the impulse method. A nice example is given in the documentation, using a limit detection function like

``````function detect_limit(u, t, integrator)
return (u[2] - 3)*(u[2] + 3)
end
``````

Should let you do a two-sided version of that example.

1 Like

So here’s an attempt at use of callbacks to change the state of the model and it is almost working.

If you modify the MWE below by changing to `MAGNITUDE = 10.0`, you’ll see a response in which the limit conditions don’t get hit. If you set `MAGNITUDE = 30.0`, the limits seem to work, but it seems to get stuck on one of them. What I expect to see is several cycles in which the position saturates at the limits, but as the driving force decays, the end looks the same as when magnitude is started smaller. I suspect this is due to some kind of deadlock condition between the callbacks, but I’m not sure how to resolve it.

BTW, I am not trying to simulate elastic collisions as in the bouncing ball demo, but an inelastic one in which if the mass is pressed against one wall or the other it will sit there until the net force is reversed as appropriate to take it off the wall (I am not conserving momentum).

``````using DifferentialEquations, Plots
plotly()

@enum state neglimit moving poslimit

NEGLIMIT = -10.0
POSLIMIT = 10.0

MAGNITUDE = 30.0

fdrive(t) = MAGNITUDE*sin(2π*0.1*t)*exp(-t/15.0)
fspring(x) = -1.0*x
fdamping(ẋ) = -.7*ẋ
force(t,x,ẋ) = fdrive(t)+fspring(x)+fdamping(ẋ)

# if we are running into the negative travel limit
neglimit_cond(u,t,integrator) = NEGLIMIT-u[2]
function neglimit_affect!(integrator)
integrator.p[1] = neglimit
integrator.u[1] = 0.0
end
nlccb = ContinuousCallback(neglimit_cond, neglimit_affect!, nothing)

# if we are running into the positive travel limit
poslimit_cond(u,t,integrator) = u[2]-POSLIMIT
function poslimit_affect!(integrator)
integrator.p[1] = poslimit
integrator.u[1] = 0.0
end
plccb = ContinuousCallback(poslimit_cond, poslimit_affect!, nothing)

# if we are trying to get off of either travel limit condition
function moving_cond(u,t,integrator)
ẋ,x = u
state, = integrator.p
state == moving   ? 0.0 : force(t,x,ẋ)
end
function moving_affect!(integrator)
integrator.p[1] = moving
end
mvccb = ContinuousCallback(moving_cond, moving_affect!)

cbs = CallbackSet(nlccb,plccb,mvccb)

function spring_mass_ODE!(du,u,p,t)
ẋ, x = u
state, mass = p
if state == moving
du[1] = 1/mass*force(t,x,ẋ)
du[2] = ẋ
else
du[1] = 0.0
du[2] = 0.0
end
end

function smsim()
tspan = [0.0, 50.0]
u0 = [0.0, 0.0]
p=[moving, 1.0 ]
prob = ODEProblem(spring_mass_ODE!, u0, tspan,p)
sol = solve(prob, Tsit5(), callback=cbs, abstol=1e-7,reltol=1e-7)
end

function plotsmsim(sol)
p1=plot(sol, vars=[1,],ylabel="velocity")
p2=plot(sol, vars=[2,],ylabel="position")
end

sol = smsim()
plotsmsim(sol)
``````

`MAGNITUDE = 10.0` response:

`MAGNITUDE = 30.0` response:

I think the problem here, if I understand your system correctly, is that you need a different condition to leave each limit stop. The force must be positive to leave the negative stop and negative to leave the positive stop. `moving_cond` just waits for force to cross 0, which is not necessarily correct, and the reason it gets stuck at the negative limit is because the spring force dominates and the force stays positive.

I tried using a `DiscreteCallback` to fix this:

``````function start_moving(u, t, integrator)
ẋ, x = u
state, = integrator.p
((state == neglimit) && (force(t, x, ẋ) > 0) ||
(state == poslimit) && (force(t, x, ẋ) < 0))
end
mvdcb = DiscreteCallback(start_moving, moving_affect!)
cbs = CallbackSet(nlccb,plccb,mvdcb)
``````

This works differently, it get stuck early but then unsticks, and I’m not sure why. I added force to the plot, it goes negative, but the systems remains stuck at the positive limit.

``````function plotsmsim(sol)
p1=plot(sol, vars=[1,],ylabel="velocity")
p2=plot(sol, vars=[2,],ylabel="position")
tspan = sol.t[1]:0.01:sol.t[end]
p3=plot(tspan, force.(tspan, sol(tspan, idxs=2), sol(tspan, idxs=1)), ylabel="force")
end
``````

Setting `dtmax=0.1` in the call to the solver does fix things, but this seems like a heavy hammer, I’m not sure how to give a gentler hint.

Thanks for the help. This is one of several things I’m working on right now and so I’m only able to get to it intermittently. I will do some troubleshooting by adding logging to the callbacks.

I now have two solutions but neither is quite right.

First, combining ContinuousCallbacks and a DiscreteCallback works but I have to set `dtmax` to a small value to get accuracy (is there an equivalent of `rootfind=true` for DiscreteCallback?)

``````using DifferentialEquations, Plots
plotly()

@enum state neglimit moving poslimit

NEGLIMIT = -10.0
POSLIMIT = 10.0

MAGNITUDE = 25.0

fdrive(t) = MAGNITUDE*sin(2π*0.1*t)*exp(-t/15.0)
fspring(x) = -1.0*x
fdamping(ẋ) = -.7*ẋ
force(t,x,ẋ) = fdrive(t)+fspring(x)+fdamping(ẋ)

# if we are running into the negative travel limit
neglimit_cond(u,t,integrator) = NEGLIMIT-u[2]
function neglimit_affect!(integrator)
integrator.p[1] = neglimit
integrator.p[2] = neglimit
integrator.u[1] = 0.0
end
nlccb = ContinuousCallback(neglimit_cond, neglimit_affect!, nothing)

# if we are running into the positive travel limit
poslimit_cond(u,t,integrator) = u[2]-POSLIMIT
function poslimit_affect!(integrator)
integrator.p[1] = poslimit
integrator.p[2] = poslimit
integrator.u[1] = 0.0
end
plccb = ContinuousCallback(poslimit_cond, poslimit_affect!, nothing)

# if we are trying to get off of neglimit travel limit condition
function moving_cond(u,t,integrator)
ẋ,x = u
state, callback_prev_state = integrator.p
if callback_prev_state == neglimit
out = force(t,x,ẋ)>0
elseif callback_prev_state == poslimit
out = force(t,x,ẋ)<0
else
out = false
end
out
end
function moving_affect!(integrator)
integrator.p[1] = moving
end
mvccb = DiscreteCallback(moving_cond, moving_affect!)

cbs = CallbackSet(nlccb,plccb,mvccb)

function spring_mass_ODE!(du,u,p,t)
ẋ, x = u
state, calback_prev_force, mass = p
if state == moving
du[1] = 1/mass*force(t,x,ẋ)
du[2] = ẋ
else
du[1] = 0.0
du[2] = 0.0
end
end

function smsim()
tspan = [0.0, 50.0]
u0 = [0.0, 0.0]
p=[moving, moving, 1.0 ]
prob = ODEProblem(spring_mass_ODE!, u0, tspan,p)
sol = solve(prob, Tsit5(), callback=cbs, abstol=1e-7,reltol=1e-7,dtmax=0.1)
end

function plotsmsim(sol)
p1=plot(sol, vars=[1,],ylabel="velocity")
p2=plot(sol, vars=[2,],ylabel="position")
end

function simplot()
sol = smsim()
plotsmsim(sol)
end
``````

The second uses two DiscreteCallbacks but has an issue in which the solutions “buzzes” when it’s against the limit stops.

``````using DifferentialEquations, Plots
plotly()

NEGLIMIT = -10.0
POSLIMIT = 10.0

MAGNITUDE = 25.0

fdrive(t) = MAGNITUDE*sin(2π*0.1*t)*exp(-t/15.0)
fspring(x) = -1.0*x
fdamping(ẋ) = -.7*ẋ
force(t,x,ẋ) = fdrive(t)+fspring(x)+fdamping(ẋ)

# if we are running into the negative travel limit
neglimit_cond(u,t,integrator) = NEGLIMIT-u[2]≥0
function neglimit_affect!(integrator)
ẋ, x = integrator.u
integrator.u[1] = 0.0
integrator.u[2] = NEGLIMIT
end
nlccb = DiscreteCallback(neglimit_cond, neglimit_affect!)

# if we are running into the positive travel limit
poslimit_cond(u,t,integrator) = u[2]-POSLIMIT≥0
function poslimit_affect!(integrator)
ẋ, x = integrator.u
integrator.u[1] = 0.0
integrator.u[2] = POSLIMIT
end
plccb = DiscreteCallback(poslimit_cond, poslimit_affect!)

cbs = CallbackSet(nlccb,plccb)

function spring_mass_ODE!(du,u,p,t)
ẋ, x = u
mass, = p
du[1] = 1/mass*force(t,x,ẋ)
du[2] = ẋ
end

function smsim()
tspan = [0.0, 50.0]
u0 = [0.0, 0.0]
p=[1.0]
prob = ODEProblem(spring_mass_ODE!, u0, tspan,p)
sol = solve(prob, Tsit5(), callback=cbs, abstol=1e-9,reltol=1e-9)
end

function plotsmsim(sol)
p1=plot(sol, vars=[1,],ylabel="velocity")
p2=plot(sol, vars=[2,],ylabel="position")
end

function simplot()
sol = smsim()
plotsmsim(sol)
end
``````

Attaching plots:
From the first simulation (with the buzzing):

and from the second:
.

It looks like that should be directional ContinuousCallbacks?

I tried that but it didn’t catch all cases.

When there is a collision, there are two possible scenarios:

• the driving force can still be in the direction of the wall and can hold it there until the driving force changes sign (the first three collisions visible in the plots), or
• the driving force was already of opposite sign at time of collision (fourth collision in plot), so the mass coasts into the wall and should immediately bounce off, with exactly one point at the limit position.

I first wrote similar code using ContinuousCallbacks, but they would not find the second scenario as there was no sign change (at least after the last callback event?) to be found.

[EDIT - writing this reply gave me another idea for how to attack it. Will post again in a while.]