# DifferentialEquations; how to simulate limit stops

**URL:** <https://discourse.julialang.org/t/differentialequations-how-to-simulate-limit-stops/44376>\
**Category:** Modelling & Simulations\
**Created:** [August 5, 2020, 11:31pm UTC](https://discourse.julialang.org/t/differentialequations-how-to-simulate-limit-stops/44376 "2020-08-05T23:31:52Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [August 5, 2020, 11:31pm UTC](https://discourse.julialang.org/t/differentialequations-how-to-simulate-limit-stops/44376/1 "2020-08-05T23:31:52Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [August 5, 2020, 11:53pm UTC](https://discourse.julialang.org/t/differentialequations-how-to-simulate-limit-stops/44376/2 "2020-08-05T23:53:03Z")

</div>

Adding a MWE:

```julia
using DifferentialEquations, Plots
plotly()

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

function spring_mass_ODE!(du,u,p,t)
	ẋ, x = u
	mass, k, kd = p
	du[1] = 1/mass*(-k*x-kd*ẋ+fdrive(t))
	du[2] = ẋ
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")
	plot(p1,p2,layout=(2,1),link=:x,leg=false)
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 6, 2020, 2:04am UTC](https://discourse.julialang.org/t/differentialequations-how-to-simulate-limit-stops/44376/3 "2020-08-06T02:04:37Z")

</div>

> [@klaff](#):
>
> 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).

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.

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [August 6, 2020, 5:06pm UTC](https://discourse.julialang.org/t/differentialequations-how-to-simulate-limit-stops/44376/4 "2020-08-06T17:06:27Z")

</div>

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

```julia
# 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 = 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*ẋ+negative_stop_force+positive_stop_force+fdrive(t))
    du[2] = ẋ
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](https://en.wikipedia.org/wiki/Collision_response#Impulse-based_contact_model) method. A nice [example](https://diffeq.sciml.ai/v6.8/features/callback_functions/#Example-1:-Bouncing-Ball-1) is given in the documentation, using a limit detection function like

```julia
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.

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [August 7, 2020, 9:30pm UTC](https://discourse.julialang.org/t/differentialequations-how-to-simulate-limit-stops/44376/5 "2020-08-07T21:30:40Z")

</div>

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

```julia
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(ẋ) = -.7*ẋ
force(t,x,ẋ) = fdrive(t)+fspring(x)+fdamping(ẋ)

# 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 = u
	state, = integrator.p
	state == moving ? 0.0 : force(t,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 = u
	state, mass = p
	if state == moving
		du[1] = 1/mass*force(t,x,ẋ)
		du[2] = ẋ
	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")
	plot(p1,p2,layout=(2,1),link=:x,leg=false)
end
	
sol = smsim()
plotsmsim(sol)

```

`MAGNITUDE = 10.0` response:  
 ![mag10](https://global.discourse-cdn.com/julialang/original/3X/d/1/d1f99e9c080b408caf4f2f3599dd62ed2db738cb.png)

`MAGNITUDE = 30.0` response:  
 ![mag30](https://global.discourse-cdn.com/julialang/original/3X/d/3/d3d2c3faa775e5acefc5980c31f126ebe697163c.png)

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [August 8, 2020, 7:51pm UTC](https://discourse.julialang.org/t/differentialequations-how-to-simulate-limit-stops/44376/6 "2020-08-08T19:51:27Z")

</div>

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:

```julia
function start_moving(u, t, integrator)                                          
    ẋ, x = u                                                                     
    state, = integrator.p                                                        
    ((state == neglimit) && (force(t, x, ẋ) > 0) ||                              
     (state == poslimit) && (force(t, 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.

```julia
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")
    plot(p1,p2,p3,layout=(3,1),link=:x,leg=false)
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.

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [August 10, 2020, 6:20pm UTC](https://discourse.julialang.org/t/differentialequations-how-to-simulate-limit-stops/44376/7 "2020-08-10T18:20:43Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [August 10, 2020, 9:26pm UTC](https://discourse.julialang.org/t/differentialequations-how-to-simulate-limit-stops/44376/8 "2020-08-10T21:26:09Z")

</div>

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

```julia
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(ẋ) = -.7*ẋ
force(t,x,ẋ) = fdrive(t)+fspring(x)+fdamping(ẋ)

# 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 = u
    state, callback_prev_state = integrator.p
    if callback_prev_state == neglimit
        out = force(t,x,ẋ)>0
    elseif callback_prev_state == poslimit
        out = force(t,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 = u
    state, calback_prev_force, mass = p
    if state == moving
        du[1] = 1/mass*force(t,x,ẋ)
        du[2] = ẋ
    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")
    plot(p1,p2,layout=(2,1),link=:x,leg=false)
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.

```julia
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(ẋ) = -.7*ẋ
force(t,x,ẋ) = fdrive(t)+fspring(x)+fdamping(ẋ)

# if we are running into the negative travel limit
neglimit_cond(u,t,integrator) = NEGLIMIT-u[2]≥0
function neglimit_affect!(integrator)
    ẋ, 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 = 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 = u
    mass, = p
    du[1] = 1/mass*force(t,x,ẋ)
    du[2] = ẋ
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")
    plot(p1,p2,layout=(2,1),link=:x,leg=false)
end

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

```

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [August 10, 2020, 9:28pm UTC](https://discourse.julialang.org/t/differentialequations-how-to-simulate-limit-stops/44376/9 "2020-08-10T21:28:17Z")

</div>

Attaching plots:  
From the first simulation (with the buzzing):  
 ![sol1](https://global.discourse-cdn.com/julialang/original/3X/c/1/c1697a43bf31a6e8427788a649da65e1bdc8fd49.png)  
and from the second:  
 ![sol2](https://global.discourse-cdn.com/julialang/original/3X/f/5/f5a8e78f86b1eb001e0271ca5adab7d92a39f59b.png).

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 10, 2020, 9:32pm UTC](https://discourse.julialang.org/t/differentialequations-how-to-simulate-limit-stops/44376/10 "2020-08-10T21:32:31Z")

</div>

> [@klaff](#):
>
> ```julia
> function moving_cond(u,t,integrator)
> ẋ,x = u
> state, callback_prev_state = integrator.p
> if callback_prev_state == neglimit
> out = force(t,x,ẋ)>0
> elseif callback_prev_state == poslimit
> out = force(t,x,ẋ)<0
> else
> out = false
> end
> out
> end
> 
> ```

It looks like that should be directional ContinuousCallbacks?

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [August 10, 2020, 10:22pm UTC](https://discourse.julialang.org/t/differentialequations-how-to-simulate-limit-stops/44376/11 "2020-08-10T22:22:17Z")

</div>

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