Coulomb and static friction with ModelingToolkit.jl and DifferentialEquations.jl

Specific Domains Modelling & Simulations
,
1

I am trying to create a model of a damped oscillator with friction. It works fine with only Coulomb friction, but all my attempts to add a static friction term have failed so far. Below is the most promising, but I can’t figure out why it fails.

The failure itself would be good to understand, but I would also like to understant the expresion in the traceback. It looks like my expression, but has an implementation of sign(x(t)) in it. I would love to be able to use that, but when I add anything with a Boolean to my expression, I get ERROR: TypeError: non-boolean (Term{Real}) used in boolean context

What is going on here? What else should I try? I have seen this problem solved by adding constraints that lock motion until the force is large enough to overcome static friction when the velocity falls to zero, Is there a way to do that?

using ModelingToolkit, DifferentialEquations

@parameters t k M σs σd vmin
@variables x(t) v(t)
@derivatives D'~t

tsgn(x) = tanh(x*10)
Ht(x) = (1 + tsgn(x)) / 2

eqs = [
    D(x) ~ v,
    D(v) ~ (Ht(abs(v) - vmin) * (k * x + σd * tsgn(v)) +
            Ht(vmin - abs(v)) * Ht(k * abs(x) - σs) * (k * x - tsgn(x) * σs)) / M
]

sys = ODESystem(eqs, t)

u0 = [x=>0, v=>10]

p = [k=>4.0, M=>0.25, σd=>0.1, σs=>1.0, vmin=>0.01]

tspan = (0.0, 50.0)

prob = ODEProblem(sys, u0, tspan, p, jac=true)

sol = solve(prob, Rodas4P());

plot(sol, vars=(0,2))

This results in the following stack trace when creating the ODEProblem.

ERROR: Failed to apply rule (*)(~(~x)) => if length(~(~x)) > N
        (*)((*)((~(~x))[1:N]...) * (*)((~(~x))[N + 1:end]...))
    else
        nothing
    end on expression 2.5 * (1 + tanh(10 * ((-1 * abs(v⦗t⦘)) + vmin))) * ((k * x⦗t⦘) + (-1 * tanh(10 * x⦗t⦘) * σs)) * (1 + (-1 * (tanh(10 * ((k * abs(x⦗t⦘)) + (-1 * σs))) ^ 2))) * ifelse(signbit(x⦗t⦘), -1, 1) * k
Stacktrace:
 [1] (::SymbolicUtils.Rule{Term{Any},SymbolicUtils.var"#term_matcher#46"{Tuple{SymbolicUtils.var"#literal_matcher#42"{typeof(*)},SymbolicUtils.var"#segment_matcher#44"{SymbolicUtils.Segment{typeof(SymbolicUtils.alwaystrue)}}}},ModelingToolkit.var"#453#455"{Int64}})(::Term{Number}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/rule.jl:142
 [2] (::SymbolicUtils.Rewriters.Chain)(::Term{Number}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/rewriters.jl:59
 [3] (::SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false})(::Term{Number}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/rewriters.jl:140
 [4] (::SymbolicUtils.Rewriters.PassThrough{SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false}})(::Term{Number}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/rewriters.jl:128
 [5] iterate at ./generator.jl:47 [inlined]
 [6] collect_to!(::Array{Term{Number},1}, ::Base.Generator{Array{Term{Number},1},SymbolicUtils.Rewriters.PassThrough{SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false}}}, ::Int64, ::Int64) at ./array.jl:732
 [7] collect_to_with_first!(::Array{Term{Number},1}, ::Term{Number}, ::Base.Generator{Array{Term{Number},1},SymbolicUtils.Rewriters.PassThrough{SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false}}}, ::Int64) at ./array.jl:710
 [8] _collect(::Array{Term{Number},1}, ::Base.Generator{Array{Term{Number},1},SymbolicUtils.Rewriters.PassThrough{SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false}}}, ::Base.EltypeUnknown, ::Base.HasShape{1}) at ./array.jl:704
 [9] collect_similar at ./array.jl:628 [inlined]
 [10] map(::SymbolicUtils.Rewriters.PassThrough{SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false}}, ::Array{Term{Number},1}) at ./abstractarray.jl:2162
 [11] (::SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false})(::Term{Number}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/rewriters.jl:138
 [12] (::SymbolicUtils.Rewriters.PassThrough{SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false}})(::Term{Number}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/rewriters.jl:128
 [13] iterate at ./generator.jl:47 [inlined]
 [14] _collect(::Array{Term,1}, ::Base.Generator{Array{Term,1},SymbolicUtils.Rewriters.PassThrough{SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false}}}, ::Base.EltypeUnknown, ::Base.HasShape{1}) at ./array.jl:699
 [15] collect_similar at ./array.jl:628 [inlined]
 [16] map(::SymbolicUtils.Rewriters.PassThrough{SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false}}, ::Array{Term,1}) at ./abstractarray.jl:2162
 [17] (::SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false})(::Term{Number}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/rewriters.jl:138
 [18] (::SymbolicUtils.Rewriters.Fixpoint{SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false}})(::Term{Number}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/rewriters.jl:101
 [19] unflatten_long_ops(::Num, ::Int64) at /home/russel/.julia/packages/ModelingToolkit/iu31N/src/build_function.jl:97
 [20] unflatten_long_ops at /home/russel/.julia/packages/ModelingToolkit/iu31N/src/build_function.jl:91 [inlined]
 [21] #62 at ./operators.jl:875 [inlined]
 [22] iterate at ./generator.jl:47 [inlined]
 [23] collect_to!(::Array{Int64,2}, ::Base.Generator{Array{Num,2},Base.var"#62#63"{typeof(ModelingToolkit.unflatten_long_ops),ModelingToolkit.var"#476#495"{Dict{Any,Any}}}}, ::Int64, ::Int64) at ./array.jl:732
 [24] collect_to_with_first!(::Array{Int64,2}, ::Int64, ::Base.Generator{Array{Num,2},Base.var"#62#63"{typeof(ModelingToolkit.unflatten_long_ops),ModelingToolkit.var"#476#495"{Dict{Any,Any}}}}, ::Int64) at ./array.jl:710
 [25] collect(::Base.Generator{Array{Num,2},Base.var"#62#63"{typeof(ModelingToolkit.unflatten_long_ops),ModelingToolkit.var"#476#495"{Dict{Any,Any}}}}) at ./array.jl:691
 [26] _build_function(::ModelingToolkit.JuliaTarget, ::Array{Num,2}, ::Array{Term{Real},1}, ::Vararg{Any,N} where N; conv::ModelingToolkit.ODEToExpr, expression::Type{T} where T, checkbounds::Bool, linenumbers::Bool, multithread::Nothing, headerfun::typeof(ModelingToolkit.addheader), outputidxs::Nothing, convert_oop::Bool, force_SA::Bool, skipzeros::Bool, fillzeros::Bool, parallel::ModelingToolkit.SerialForm, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/russel/.julia/packages/ModelingToolkit/iu31N/src/build_function.jl:254
 [27] build_function(::Array{Num,2}, ::Vararg{Any,N} where N; target::ModelingToolkit.JuliaTarget, kwargs::Base.Iterators.Pairs{Symbol,Any,NTuple{5,Symbol},NamedTuple{(:conv, :expression, :checkbounds, :linenumbers, :parallel),Tuple{ModelingToolkit.ODEToExpr,DataType,Bool,Bool,ModelingToolkit.SerialForm}}}) at /home/russel/.julia/packages/ModelingToolkit/iu31N/src/build_function.jl:48
 [28] generate_jacobian(::ODESystem, ::Array{Term{Real},1}, ::Array{Sym{ModelingToolkit.Parameter{Real}},1}; simplify::Bool, sparse::Bool, kwargs::Base.Iterators.Pairs{Symbol,Any,NTuple{4,Symbol},NamedTuple{(:expression, :checkbounds, :linenumbers, :parallel),Tuple{DataType,Bool,Bool,ModelingToolkit.SerialForm}}}) at /home/russel/.julia/packages/ModelingToolkit/iu31N/src/systems/diffeqs/abstractodesystem.jl:67
 [29] ODEFunction{true,F,TMM,Ta,Tt,TJ,JVP,VJP,JP,SP,TW,TWt,TPJ,S,TCV} where TCV where S where TPJ where TWt where TW where SP where JP where VJP where JVP where TJ where Tt where Ta where TMM where F(::ODESystem, ::Array{Term{Real},1}, ::Array{Sym{ModelingToolkit.Parameter{Real}},1}, ::Array{Int64,1}; version::Nothing, tgrad::Bool, jac::Bool, eval_expression::Bool, sparse::Bool, simplify::Bool, kwargs::Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol,Symbol},NamedTuple{(:checkbounds, :linenumbers, :parallel),Tuple{Bool,Bool,ModelingToolkit.SerialForm}}}) at /home/russel/.julia/packages/ModelingToolkit/iu31N/src/systems/diffeqs/abstractodesystem.jl:147
 [30] ODEProblem{true,tType,isinplace,P,F,K,PT} where PT where K where F where P where isinplace where tType(::ODESystem, ::Array{Pair{Num,Int64},1}, ::Tuple{Float64,Float64}, ::Array{Pair{Num,Float64},1}; version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::ModelingToolkit.SerialForm, eval_expression::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/russel/.julia/packages/ModelingToolkit/iu31N/src/systems/diffeqs/abstractodesystem.jl:277
 [31] #ODEProblem#206 at /home/russel/.julia/packages/ModelingToolkit/iu31N/src/systems/diffeqs/abstractodesystem.jl:240 [inlined]
 [32] top-level scope at REPL[17]:1
caused by [exception 1]
MethodError: no method matching *(::Term{Real}, ::Symbol)
Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...) at operators.jl:538
  *(::LightGraphs.LinAlg.Noop, ::Any) at /home/russel/.julia/packages/LightGraphs/HsNig/src/linalg/graphmatrices.jl:226
  *(::ChainRulesCore.One, ::Any) at /home/russel/.julia/packages/ChainRulesCore/UayCG/src/differential_arithmetic.jl:98
  ...
Stacktrace:
 [1] (::ModelingToolkit.var"#453#455"{Int64})(::Base.ImmutableDict{Symbol,Any}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/rule.jl:274
 [2] macro expansion at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/utils.jl:11 [inlined]
 [3] (::SymbolicUtils.var"#34#35"{ModelingToolkit.var"#453#455"{Int64}})(::Base.ImmutableDict{Symbol,Any}, ::Int64) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/rule.jl:139
 [4] (::SymbolicUtils.var"#loop#47"{SymbolicUtils.var"#34#35"{ModelingToolkit.var"#453#455"{Int64}}})(::SymbolicUtils.LL{Array{Any,1}}, ::Base.ImmutableDict{Symbol,Any}, ::SymbolicUtils.LL{Tuple{SymbolicUtils.var"#literal_matcher#42"{typeof(*)},SymbolicUtils.var"#segment_matcher#44"{SymbolicUtils.Segment{typeof(SymbolicUtils.alwaystrue)}}}}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/matchers.jl:96
 [5] (::SymbolicUtils.var"#45#48"{SymbolicUtils.LL{Array{Any,1}},SymbolicUtils.LL{Tuple{SymbolicUtils.var"#literal_matcher#42"{typeof(*)},SymbolicUtils.var"#segment_matcher#44"{SymbolicUtils.Segment{typeof(SymbolicUtils.alwaystrue)}}}}})(::Base.ImmutableDict{Symbol,Any}, ::Int64) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/matchers.jl:101
 [6] (::SymbolicUtils.var"#segment_matcher#44"{SymbolicUtils.Segment{typeof(SymbolicUtils.alwaystrue)}})(::SymbolicUtils.var"#45#48"{SymbolicUtils.LL{Array{Any,1}},SymbolicUtils.LL{Tuple{SymbolicUtils.var"#literal_matcher#42"{typeof(*)},SymbolicUtils.var"#segment_matcher#44"{SymbolicUtils.Segment{typeof(SymbolicUtils.alwaystrue)}}}}}, ::SymbolicUtils.LL{Array{Any,1}}, ::Base.ImmutableDict{Symbol,Any}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/matchers.jl:74
 [7] (::SymbolicUtils.var"#loop#47"{SymbolicUtils.var"#34#35"{ModelingToolkit.var"#453#455"{Int64}}})(::SymbolicUtils.LL{Array{Any,1}}, ::Base.ImmutableDict{Symbol,Any}, ::SymbolicUtils.LL{Tuple{SymbolicUtils.var"#literal_matcher#42"{typeof(*)},SymbolicUtils.var"#segment_matcher#44"{SymbolicUtils.Segment{typeof(SymbolicUtils.alwaystrue)}}}}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/matchers.jl:100
 [8] (::SymbolicUtils.var"#45#48"{Term{Number},Tuple{SymbolicUtils.var"#literal_matcher#42"{typeof(*)},SymbolicUtils.var"#segment_matcher#44"{SymbolicUtils.Segment{typeof(SymbolicUtils.alwaystrue)}}}})(::Base.ImmutableDict{Symbol,Any}, ::Int64) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/matchers.jl:101
 [9] (::SymbolicUtils.var"#literal_matcher#42"{typeof(*)})(::SymbolicUtils.var"#45#48"{Term{Number},Tuple{SymbolicUtils.var"#literal_matcher#42"{typeof(*)},SymbolicUtils.var"#segment_matcher#44"{SymbolicUtils.Segment{typeof(SymbolicUtils.alwaystrue)}}}}, ::Term{Number}, ::Base.ImmutableDict{Symbol,Any}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/matchers.jl:10
 [10] (::SymbolicUtils.var"#loop#47"{SymbolicUtils.var"#34#35"{ModelingToolkit.var"#453#455"{Int64}}})(::Term{Number}, ::Base.ImmutableDict{Symbol,Any}, ::Tuple{SymbolicUtils.var"#literal_matcher#42"{typeof(*)},SymbolicUtils.var"#segment_matcher#44"{SymbolicUtils.Segment{typeof(SymbolicUtils.alwaystrue)}}}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/matchers.jl:100
 [11] (::SymbolicUtils.var"#term_matcher#46"{Tuple{SymbolicUtils.var"#literal_matcher#42"{typeof(*)},SymbolicUtils.var"#segment_matcher#44"{SymbolicUtils.Segment{typeof(SymbolicUtils.alwaystrue)}}}})(::SymbolicUtils.var"#34#35"{ModelingToolkit.var"#453#455"{Int64}}, ::Tuple{Term{Number}}, ::Base.ImmutableDict{Symbol,Any}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/matchers.jl:105
 [12] (::SymbolicUtils.Rule{Term{Any},SymbolicUtils.var"#term_matcher#46"{Tuple{SymbolicUtils.var"#literal_matcher#42"{typeof(*)},SymbolicUtils.var"#segment_matcher#44"{SymbolicUtils.Segment{typeof(SymbolicUtils.alwaystrue)}}}},ModelingToolkit.var"#453#455"{Int64}})(::Term{Number}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/rule.jl:137
 [13] (::SymbolicUtils.Rewriters.Chain)(::Term{Number}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/rewriters.jl:59
 [14] (::SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false})(::Term{Number}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/rewriters.jl:140
 [15] (::SymbolicUtils.Rewriters.PassThrough{SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false}})(::Term{Number}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/rewriters.jl:128
 [16] iterate at ./generator.jl:47 [inlined]
 [17] collect_to!(::Array{Term{Number},1}, ::Base.Generator{Array{Term{Number},1},SymbolicUtils.Rewriters.PassThrough{SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false}}}, ::Int64, ::Int64) at ./array.jl:732
 [18] collect_to_with_first!(::Array{Term{Number},1}, ::Term{Number}, ::Base.Generator{Array{Term{Number},1},SymbolicUtils.Rewriters.PassThrough{SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false}}}, ::Int64) at ./array.jl:710
 [19] _collect(::Array{Term{Number},1}, ::Base.Generator{Array{Term{Number},1},SymbolicUtils.Rewriters.PassThrough{SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false}}}, ::Base.EltypeUnknown, ::Base.HasShape{1}) at ./array.jl:704
 [20] collect_similar at ./array.jl:628 [inlined]
 [21] map(::SymbolicUtils.Rewriters.PassThrough{SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false}}, ::Array{Term{Number},1}) at ./abstractarray.jl:2162
 [22] (::SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false})(::Term{Number}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/rewriters.jl:138
 [23] (::SymbolicUtils.Rewriters.PassThrough{SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false}})(::Term{Number}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/rewriters.jl:128
 [24] iterate at ./generator.jl:47 [inlined]
 [25] _collect(::Array{Term,1}, ::Base.Generator{Array{Term,1},SymbolicUtils.Rewriters.PassThrough{SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false}}}, ::Base.EltypeUnknown, ::Base.HasShape{1}) at ./array.jl:699
 [26] collect_similar at ./array.jl:628 [inlined]
 [27] map(::SymbolicUtils.Rewriters.PassThrough{SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false}}, ::Array{Term,1}) at ./abstractarray.jl:2162
 [28] (::SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false})(::Term{Number}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/rewriters.jl:138
 [29] (::SymbolicUtils.Rewriters.Fixpoint{SymbolicUtils.Rewriters.Walk{:post,SymbolicUtils.Rewriters.Chain,false}})(::Term{Number}) at /home/russel/.julia/packages/SymbolicUtils/K1sGK/src/rewriters.jl:101
 [30] unflatten_long_ops(::Num, ::Int64) at /home/russel/.julia/packages/ModelingToolkit/iu31N/src/build_function.jl:97
 [31] unflatten_long_ops at /home/russel/.julia/packages/ModelingToolkit/iu31N/src/build_function.jl:91 [inlined]
 [32] #62 at ./operators.jl:875 [inlined]
 [33] iterate at ./generator.jl:47 [inlined]
 [34] collect_to!(::Array{Int64,2}, ::Base.Generator{Array{Num,2},Base.var"#62#63"{typeof(ModelingToolkit.unflatten_long_ops),ModelingToolkit.var"#476#495"{Dict{Any,Any}}}}, ::Int64, ::Int64) at ./array.jl:732
 [35] collect_to_with_first!(::Array{Int64,2}, ::Int64, ::Base.Generator{Array{Num,2},Base.var"#62#63"{typeof(ModelingToolkit.unflatten_long_ops),ModelingToolkit.var"#476#495"{Dict{Any,Any}}}}, ::Int64) at ./array.jl:710
 [36] collect(::Base.Generator{Array{Num,2},Base.var"#62#63"{typeof(ModelingToolkit.unflatten_long_ops),ModelingToolkit.var"#476#495"{Dict{Any,Any}}}}) at ./array.jl:691
 [37] _build_function(::ModelingToolkit.JuliaTarget, ::Array{Num,2}, ::Array{Term{Real},1}, ::Vararg{Any,N} where N; conv::ModelingToolkit.ODEToExpr, expression::Type{T} where T, checkbounds::Bool, linenumbers::Bool, multithread::Nothing, headerfun::typeof(ModelingToolkit.addheader), outputidxs::Nothing, convert_oop::Bool, force_SA::Bool, skipzeros::Bool, fillzeros::Bool, parallel::ModelingToolkit.SerialForm, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/russel/.julia/packages/ModelingToolkit/iu31N/src/build_function.jl:254
 [38] build_function(::Array{Num,2}, ::Vararg{Any,N} where N; target::ModelingToolkit.JuliaTarget, kwargs::Base.Iterators.Pairs{Symbol,Any,NTuple{5,Symbol},NamedTuple{(:conv, :expression, :checkbounds, :linenumbers, :parallel),Tuple{ModelingToolkit.ODEToExpr,DataType,Bool,Bool,ModelingToolkit.SerialForm}}}) at /home/russel/.julia/packages/ModelingToolkit/iu31N/src/build_function.jl:48
 [39] generate_jacobian(::ODESystem, ::Array{Term{Real},1}, ::Array{Sym{ModelingToolkit.Parameter{Real}},1}; simplify::Bool, sparse::Bool, kwargs::Base.Iterators.Pairs{Symbol,Any,NTuple{4,Symbol},NamedTuple{(:expression, :checkbounds, :linenumbers, :parallel),Tuple{DataType,Bool,Bool,ModelingToolkit.SerialForm}}}) at /home/russel/.julia/packages/ModelingToolkit/iu31N/src/systems/diffeqs/abstractodesystem.jl:67
 [40] ODEFunction{true,F,TMM,Ta,Tt,TJ,JVP,VJP,JP,SP,TW,TWt,TPJ,S,TCV} where TCV where S where TPJ where TWt where TW where SP where JP where VJP where JVP where TJ where Tt where Ta where TMM where F(::ODESystem, ::Array{Term{Real},1}, ::Array{Sym{ModelingToolkit.Parameter{Real}},1}, ::Array{Int64,1}; version::Nothing, tgrad::Bool, jac::Bool, eval_expression::Bool, sparse::Bool, simplify::Bool, kwargs::Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol,Symbol},NamedTuple{(:checkbounds, :linenumbers, :parallel),Tuple{Bool,Bool,ModelingToolkit.SerialForm}}}) at /home/russel/.julia/packages/ModelingToolkit/iu31N/src/systems/diffeqs/abstractodesystem.jl:147
 [41] ODEProblem{true,tType,isinplace,P,F,K,PT} where PT where K where F where P where isinplace where tType(::ODESystem, ::Array{Pair{Num,Int64},1}, ::Tuple{Float64,Float64}, ::Array{Pair{Num,Float64},1}; version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::ModelingToolkit.SerialForm, eval_expression::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/russel/.julia/packages/ModelingToolkit/iu31N/src/systems/diffeqs/abstractodesystem.jl:277
 [42] #ODEProblem#206 at /home/russel/.julia/packages/ModelingToolkit/iu31N/src/systems/diffeqs/abstractodesystem.jl:240 [inlined]
 [43] top-level scope at REPL[17]:1
2

It looks like simplify failed, and @shashi would be the one to take a look at it. Thanks for the example: that should be very useful.

3

You can use IfElse.ifelse(x,y,z) as a stand-in for the standard Core.ifelse and it’ll be something that is symbolically recognized. We may need to fix a few derivative rules to do this.

4

Ah, IfElse is a great trick, I didn’t know about that one.

I almost got things working using that idea and then I realized this is not a good model of friction anyway, since it sets D(v) ~ 0, not v ~ 0, when the static friction becomes important. I need to think more about what I really mean here. Maybe I can do something fancy with callbacks.

5

OK, I got the bad model working with IfElse.


isign(x) = IfElse.ifelse(signbit(x), -1, 1)

eqs = [
    D(x) ~ v,
    D(v) ~ IfElse.ifelse(signbit(abs(v) - vmin),
                         IfElse.ifelse(signbit(k*abs(x) - σs),
                                       0,
                                       -k * x + isign(x) * σs),
                         -k * x - σd * isign(v)) / M
]

When I try to express this more naturally with <, I get

ERROR: MethodError: no method matching isless(::Term{Number}, ::Sym{Real})

SymbolicUtils only defines isless for arguments of matching type, Is there some way to turn that Term{Number}into a Term{Real}?

6

This is probably worth an issue.

7

Hmm, this already looks pretty well thought out. Number was explicitly skipped for all comparisons besides equality, and it does allow isless for T<:Real. Looking at the tests, I think this was to leave out Complex.

8

It does not define it on Number because there’s no way to compare complex numbers but maybe we should allow Number x Real and Real x Number

9

https://github.com/SciML/ModelingToolkit.jl/pull/648 this should fix the issue in the OP

10

That works, thanks! I’m surprised how similar the two solutions are, I would have expected the smooth approximation to have a much more noticeable effect.

11

The above code does not seem to be working in Win10 Julia 1.5.2. Tried to add ModelingToolkit#master to get latest fixes. I always get the following error in the eqs definition:

ERROR: MethodError: no method matching signbit(::Operation)

Thank you.

12

Operation shouldn’t exist anymore. Check ]st: do you have MTK v4? This is a very recent release and a very recent fix.

1 Like
13

@ChrisRackauckas, thank you. After a lengthy package maintenance/hygiene work, managed to update from MTK v3.1 to v3.20.1. The code runs with no errors now.

PS: How can we update to v4?

14

You probably need some updates which may not be released yet. Some downstream stuff are still in the release process (it’s that new still)

1 Like
15

Thanks for sharing this very interesting physics problem. One observation regarding the friction model equations, if you will.
One could define the static friction region as k*abs(x) < σs*N where the spring force is cancelled (N is the normal force). And the kinematic friction region over k*abs(x) >= σs*N. The equations simplify then to:

const g = 9.8;  # m/s^2
eqs = [
    D(x) ~ v,
    D(v) ~
        IfElse.ifelse(signbit.(k*abs(x) - σs*M*g), 0,  -k*x -σd*M*g*isign(v)) / M,
]

Also, the equations using the smooth step function trick would then become:

eqs1 = [
    D(x) ~ v,
    D(v) ~
        ( Ht(k*abs(x) - σs*M*g) * (-k*x - σd*M*g*tsgn(v)) ) / M,
]

Best regards.

3 Likes