# ModelingToolkit.jl and differential equations with step/pulse inputs

**URL:** <https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272>\
**Category:** Modelling & Simulations\
**Tags:** diffeq\
**Created:** [May 27, 2020, 4:27pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272 "2020-05-27T16:27:58Z")\
**Posts on this page:** 18\
**Page:** 1

<div class="post-metadata">

**Author:** ![DR59](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dr59/32/10087_2.png) [@DR59](https://discourse.julialang.org/u/DR59)\
**Post date:** [May 27, 2020, 4:27pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272/1 "2020-05-27T16:27:58Z")

</div>

In lots of differential equations applications, we want to solve a ‘forced’ DE whose evolution depends upon some input u(t). For instance we can consider a harmonic oscillator with a forcing term. The equation below makes a function create(), that takes an input u(t), and outputs an ODEProblem of the form

\dot{x}(t) = f(x,p, u(t)),  
representing a harmonic oscillator being forced by your choice of input u(t), with parameters p.

```julia
function create(input)
    
  @parameters t
  @parameters k,c,m
  @derivatives D'~t
  @variables pos(t) vel(t) inp(t)

  eqs = [D(pos) ~ vel,
        D(vel) ~ input(t) - (1/m)*(c*vel + k*pos),
  ]
  ps = [k,c,m] .=> [1.,4.,1.] 
  ics = [pos, vel] .=> [1.,0.]
  od = ODESystem(eqs, t, first.(ics),first.(ps))
  tspan = (0.,100.)
  prob = ODEProblem(od, ics, tspan, ps)
  return prob
end

```

If I take ‘nice’ inputs, everything works great. For instance,  
`prob = create(sin)`

However, in a lot of important applications (eg pharmacokinetic modelling), inputs come as nondifferentiable functions, such as pulses or steps. E.g.

`nasty(t) = t>1 ? 1 : 0`

Now I can understand that functions with if statements and the like, are not convertable to symbolic code by ModelingToolkit.jl. Indeed, the code:

```julia
@parameters t
nasty(t)

```

gives an error:

> ERROR: TypeError: non-boolean (Operation) used in boolean context

But it would be really nice if ModelingToolkit.jl could generate one of its’ awesome, compiled vector fields from eqs using inputs like nasty(t), without ever forming an (impossible) symbolic representation of the nasty bit. Basically I want some way of doing

`prob = create(nasty)`

Since a function like nasty(t) depends only on t, we don’t need to do anything symbolic (jacobians etc) with it when converting the overall system of equations eqs into the vector field of an ODEProblem.

Is there any way of doing something like this using ModelingToolkit.jl? Is it even reasonable?

Many thanks, and also a big thumbs up to everyone working on ModelingToolkit.jl. Although I’m not enough of a programmer to directly contribute, I’m really enjoying seeing it develop: it’s very cool. Hopefully I can contribute a few implementations of ODE models from the systems biology literature in the fullness of time.

---

<div class="post-metadata">

**Author:** ![asinghvi17](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/asinghvi17/32/8272_2.png) [@asinghvi17](https://discourse.julialang.org/u/asinghvi17)\
**Post date:** [May 28, 2020, 10:43am UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272/2 "2020-05-28T10:43:19Z")

</div>

A quick-and-dirty solution to the problem might be to overload the `ifelse` function for Operations, or create an `@ifelse` macro which does a similar thing. For now, you could use `min` and `max` to “simulate” an `if` statement (or a ternary operator)…

---

<div class="post-metadata">

**Author:** ![DR59](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dr59/32/10087_2.png) [@DR59](https://discourse.julialang.org/u/DR59)\
**Post date:** [May 28, 2020, 11:00am UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272/3 "2020-05-28T11:00:50Z")

</div>

Thanks for the reply!

So replacing my ternary operator with a minmax-style (equivalent) operator works! Pulses would get somewhat annoying with this type of syntax I imagine.

I guess I’m not sure what you mean by overloading the ifelse function. Is there a page of documentation / a file in the source code that implements something similar? Cheers!

---

<div class="post-metadata">

**Author:** ![asinghvi17](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/asinghvi17/32/8272_2.png) [@asinghvi17](https://discourse.julialang.org/u/asinghvi17)\
**Post date:** [May 28, 2020, 11:31am UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272/4 "2020-05-28T11:31:40Z")

</div>

Ah, on closer inspection it looks like someone tried that:

> <https://github.com/SciML/ModelingToolkit.jl/blob/89cf39a288f38f2e5f7b8bfbfe11a909c3651f1e/src/function_registration.jl#L59-L60>

Unfortunately, `ifelse` is not in `Base`, but `Core`. This means that it’s a builtin function and cannot be extended, so that idea is moot. _However_, it is possible to do:

```julia
julia> @inline function if_else(condition::Bool, @nospecialize(trueval), @nospecialize(falseval))
           return ifelse(condition, trueval, falseval)
       end
if_else (generic function with 1 method)

julia> ModelingToolkit.@register if_else(x, trueval, falseval)
if_else (generic function with 8 methods)

julia> if_else |> methods
# 8 methods for generic function "if_else":
[1] if_else(condition::Bool, trueval, falseval) in Main at REPL[13]:2
[2] if_else(x::Expression, trueval::Expression, falseval::Expression) in Main
[3] if_else(x::Number, trueval::Expression, falseval::Expression) in Main
[4] if_else(x::Expression, trueval::Number, falseval::Expression) in Main
[5] if_else(x::Number, trueval::Number, falseval::Expression) in Main
[6] if_else(x::Expression, trueval::Expression, falseval::Number) in Main
[7] if_else(x::Number, trueval::Expression, falseval::Number) in Main
[8] if_else(x::Expression, trueval::Number, falseval::Number) in Main

julia> nasty(i) = if_else(i > 1, 1, 0)
nasty (generic function with 1 method)

julia> @variables t
(t,)

julia> nasty(t)
if_else(t > 1, 1, 0)

```

This _should_ work in a differential equation.

Indeed,

```julia
julia> prob = create(nasty)
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 100.0)
u0: [1.0, 0.0]

```

cc @ChrisRackauckas

Edit: Oh, I forgot about this…functions are eval’ed in ModelingToolkit’s scope, so you would get

```julia
var"##MTIIPVar#258"[2] = (ModelingToolkit).:-((ModelingToolkit).if_else((ModelingToolkit).:>(t, 1), 1, 0), (ModelingToolkit).:*((ModelingToolkit).:/(1, m), (ModelingToolkit).:+((ModelingToolkit).:*(c, vel), (ModelingToolkit).:*(k, pos))))

```

this IR. Note the `ModelingToolkit.if_else` - you may have to `@eval ModelingToolkit using Main: if_else`.

---

<div class="post-metadata">

**Author:** ![jonniedie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonniedie/32/12842_2.png) [@jonniedie](https://discourse.julialang.org/u/jonniedie)\
**Post date:** [May 28, 2020, 12:42pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272/5 "2020-05-28T12:42:55Z")

</div>

Not the most general solution to your problem, but:

```julia
less_nasty(t) = 1 / (1 + exp(1e3 * (1 - t)))

```

Or generalizing:

```julia
less_nasty(t, nastiness, step_time) = 1 / (1 + exp(nastiness * (step_time - t)))

```

With a closure:

```julia
less_nasty(nastiness, step_time) = t -> less_nasty(t, nastiness, step_time)

```

---

<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:** [May 28, 2020, 1:40pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272/6 "2020-05-28T13:40:02Z")

</div>

> [@asinghvi17](#):
>
> Unfortunately, `ifelse` is not in `Base` , but `Core` . This means that it’s a builtin function and cannot be extended, so that idea is moot

Yeah… this is … something that hasn’t worked out so well. We need a better solution on this problem. Maybe just make `if_else` and go with it, but that’s not great.

---

<div class="post-metadata">

**Author:** ![DR59](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dr59/32/10087_2.png) [@DR59](https://discourse.julialang.org/u/DR59)\
**Post date:** [May 28, 2020, 6:57pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272/7 "2020-05-28T18:57:00Z")

</div>

Nice! For steps that works perfectly. I guess with a bit of care you could also do pulses in a similar manner. Useful as a backup hack.

---

<div class="post-metadata">

**Author:** ![DR59](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dr59/32/10087_2.png) [@DR59](https://discourse.julialang.org/u/DR59)\
**Post date:** [May 28, 2020, 7:21pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272/8 "2020-05-28T19:21:11Z")

</div>

Thanks! That clarifies things a lot.

Unfortunately, while

```julia
prob = create(nasty)

```

works for your construction of nasty(),

```julia
prob |> solve

```

returns a long error message that I’m not sure how to interpret. I’m happy to use continuous hacks for now. Thanks again!

---

<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:** [May 28, 2020, 7:27pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272/9 "2020-05-28T19:27:12Z")

</div>

> [@DR59](#):
>
> returns a long error message that I’m not sure how to interpret.

post it?

> [@DR59](#):
>
> Nice! For steps that works perfectly. I guess with a bit of care you could also do pulses in a similar manner. Useful as a backup hack.

It’s probably even better, because it’ll jumps actually cause order loss in ODE solvers, so relaxing them will improve simulations.

---

<div class="post-metadata">

**Author:** ![DR59](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dr59/32/10087_2.png) [@DR59](https://discourse.julialang.org/u/DR59)\
**Post date:** [May 28, 2020, 8:02pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272/10 "2020-05-28T20:02:28Z")

</div>

By long I mean really long. Doesn’t fit on my terminal long! I’ll post two snippets? The first is below, and looks like a portion of an assembly language representation of my function. Most of the error message looks similar. The second is the end of the error message, which might be more enlightening?

Is that enough?  
Cheers

```julia
 begin
            $(Expr(:inbounds, true))
            local var"#1060#val" = begin
                        #= /.julia/packages/ModelingToolkit/ErYw5/src/build_function.jl:252 =#
                        let (pos, vel, k, c, m, t) = (var"##MTKArg#412"[1], var"##MTKArg#412"[2], var"##MTKArg#413"[1], var"##MTKArg#413"[2], var"##MTKArg#413"[3], var"##MTKArg#414")
                            begin
                                var"##MTIIPVar#416"[1] = vel
                                var"##MTIIPVar#416"[2] = (ModelingToolkit).:-((ModelingToolkit).if_else((ModelingToolkit).:>(t, 1), 1, 0), (ModelingToolkit).:*((ModelingToolkit).:/(1, m), (ModelingToolkit).:+((ModelingToolkit).:*(c, vel), (ModelingToolkit).:*(k, pos))))
                            end
                        end
                    end
            $(Expr(:inbounds, :((ModelingToolkit).pop)))
            var"#1060#val"
        end
        #= /.julia/packages/ModelingToolkit/ErYw5/src/build_function.jl:69 =#
        (ModelingToolkit).nothing
    end
end),:function}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Array{Symbol,1},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,ModelingToolkit.var"#f#86"{GeneralizedGenerated.NGG.RuntimeFn{TypeEncoding(list(##MTKArg#412, ##MTKArg#413, ##MTKArg#414)),TypeEncoding(nil(GeneralizedGenerated.NGG.Argument)),TypeEncoding(begin

```

Second

```julia
 end),:function}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Array{Symbol,1},Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}) at /.julia/packages/DiffEqBase/TiYWr/src/solve.jl:66
 [18] top-level scope at REPL[30]:1
 [19] eval(::Module, ::Any) at ./boot.jl:331
 [20] eval_user_input(::Any, ::REPL.REPLBackend) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
 [21] run_backend(::REPL.REPLBackend) at /.julia/packages/Revise/MgvIv/src/Revise.jl:1023
 [22] top-level scope at none:0

```

---

<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:** [May 28, 2020, 8:08pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272/11 "2020-05-28T20:08:59Z")

</div>

That’s not the right spot.

---

<div class="post-metadata">

**Author:** ![DR59](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dr59/32/10087_2.png) [@DR59](https://discourse.julialang.org/u/DR59)\
**Post date:** [May 28, 2020, 10:13pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272/12 "2020-05-28T22:13:12Z")

</div>

I think the right spot is the beginning of the error message. But I can’t see it because the output saturates my terminal. Even the portion of the error message displayed in my terminal can’t be shown here, as it breaks the character limit. I could send a .txt file to you with the message? Or just the script that generated it, which is below. Sorry about that.

```julia
using ModelingToolkit
using DifferentialEquations

@inline function if_else(condition::Bool, @nospecialize(trueval), @nospecialize(falseval))
           return ifelse(condition, trueval, falseval)
end

ModelingToolkit.@register if_else(x, trueval, falseval)

nasty(i) = if_else(i > 1, 1, 0)
@variables t
nasty(t)

function create(input)
    
  @parameters t
  @parameters k,c,m
  @derivatives D'~t
  @variables pos(t) vel(t) inp(t)

  eqs = [D(pos) ~ vel,
        D(vel) ~ input(t) - (1/m)*(c*vel + k*pos),
  ]

  ps = [k,c,m] .=> [1.,4.,1.] 
  ics = [pos, vel] .=> [1.,0.]
  od = ODESystem(eqs, t, first.(ics),first.(ps))
  tspan = (0.,100.)
  prob = ODEProblem(od, ics, tspan, ps)
  return prob
end

prob = create(nasty)
sol = solve(prob)

```

---

<div class="post-metadata">

**Author:** ![asinghvi17](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/asinghvi17/32/8272_2.png) [@asinghvi17](https://discourse.julialang.org/u/asinghvi17)\
**Post date:** [May 28, 2020, 11:25pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272/13 "2020-05-28T23:25:44Z")

</div>

For me, this works:

```julia
using ModelingToolkit
using DifferentialEquations

@inline function if_else(condition::Bool, @nospecialize(trueval), @nospecialize(falseval))
           return ifelse(condition, trueval, falseval)
end

ModelingToolkit.@register if_else(x, trueval, falseval)
# This is needed to bring if_else into MTK's scope.
@eval ModelingToolkit using ..Main: if_else

nasty(i) = if_else(i > 1, 1, 0)
@variables t
nasty(t)

function create(input)
    
  @parameters t
  @parameters k,c,m
  @derivatives D'~t
  @variables pos(t) vel(t) inp(t)

  eqs = [D(pos) ~ vel,
        D(vel) ~ input(t) - (1/m)*(c*vel + k*pos),
  ]

  ps = [k,c,m] .=> [1.,4.,1.] 
  ics = [pos, vel] .=> [1.,0.]
  od = ODESystem(eqs, t, first.(ics),first.(ps))
  tspan = (0.,100.)
  prob = ODEProblem(od, ics, tspan, ps)
  return prob
end

prob = create(nasty)
sol = solve(prob)

```

![sol](https://global.discourse-cdn.com/julialang/original/3X/6/a/6af2a54f0eef58eea42e00db7d5b5106bc3a92e3.png)

---

<div class="post-metadata">

**Author:** ![DR59](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dr59/32/10087_2.png) [@DR59](https://discourse.julialang.org/u/DR59)\
**Post date:** [May 29, 2020, 9:34pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272/14 "2020-05-29T21:34:21Z")

</div>

Ignore me, sorry. It works fine! My REPL was dirty I think. When I refreshed it, there was no error 🙂

---

<div class="post-metadata">

**Author:** ![DR59](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dr59/32/10087_2.png) [@DR59](https://discourse.julialang.org/u/DR59)\
**Post date:** [May 29, 2020, 9:35pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272/15 "2020-05-29T21:35:00Z")

</div>

Indeed, it works for me too. It was a stupid bug due to (I think) a dirty REPL. Thanks a lot for all your help!

---

<div class="post-metadata">

**Author:** ![PikaPei](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pikapei/32/16872_2.png) [@PikaPei](https://discourse.julialang.org/u/PikaPei)\
**Post date:** [August 5, 2020, 2:59am UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272/16 "2020-08-05T02:59:45Z")

</div>

Hello,

I just wonder is there a better solution now?

I have a similar problem like the original post, but I need a square wave input function.

I have tried the method discussed above.

```julia
@inline function if_else(condition::Bool, @nospecialize(trueval), @nospecialize(falseval))
           return ifelse(condition, trueval, falseval)
end

ModelingToolkit.@register if_else(x, trueval, falseval)
@eval ModelingToolkit using ..Main: if_else

```

If the input is a step function, it works well. Like this.

```julia
nasty(i) = if_else(i > 1, 1, 0)

```

However, if I limit the range.

```julia
nasty(i) = if_else(1 < i < 2, 1, 0)

```

It throws an error.

```julia
TypeError: non-boolean (Operation) used in boolean context

```

I’m excited to use ModelingToolkit.jl to do some simulations, but I don’t know how to solve this. If someone would help, I’d be very grateful!

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [August 8, 2020, 6:16pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272/17 "2020-08-08T18:16:35Z")

</div>

> [@PikaPei](#):
>
> `nasty(i) = if_else(i > 1, 1, 0)`

What about if you define

```julia
step(x) = if_else( a < x, 1, 0) * if_else(x < b, 1, 0)

```

---

<div class="post-metadata">

**Author:** ![PikaPei](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pikapei/32/16872_2.png) [@PikaPei](https://discourse.julialang.org/u/PikaPei)\
**Post date:** [August 9, 2020, 4:42pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-and-differential-equations-with-step-pulse-inputs/40272/18 "2020-08-09T16:42:52Z")

</div>

Yes, it works!

It is really a simple solution, many thanks!