State selection

[YingboMa]

Without state selection, ODAEProblem w/ Tsit5 gives terrible numerical drift on the pendulum model.

using ModelingToolkit, LinearAlgebra, OrdinaryDiffEq

function Pendulum(;name, g=9.81, r=1.)
    gval, rval = g, r
    @parameters g r
    @variables x(t) y(t) T(t)
    D2 = Differential(t)^2
    eqs = [
           D2(x) ~ T * x
           D2(y) ~ T * y - g
           x^2 + y^2 ~ r
          ]

    ODESystem(eqs, t, [x, y, T], [g, r]; name=name, defaults=(g=>gval, r=>rval, x=>1, y=>0, D(x)=>0, D(y)=>0.0, T=>0))
end

@parameters t
D = Differential(t)
@named sys = Pendulum(;g=9.81, r=1)
pendulum_sys = structural_simplify(ode_order_lowering(flatten(sys)))
prob = ODAEProblem(pendulum_sys, [], (0, 100.0))
sol = solve(prob, Tsit5())
x = @nonamespace sys.x
y = @nonamespace sys.y
@. sol[x]^2 + sol[y]^2

and

julia> using Plots

julia> plot(sol, vars=(x, y), xlab="x", ylab="y", dpi=300, title="Numerical drift")

gives

[anandijain]

could you post a reference? i still dont know what state selection is lol

[kabdelhak]

could you post a reference? i still dont know what state selection is lol

https://www.researchgate.net/publication/235324214_Index_Reduction_in_Differential-Algebraic_Equations_Using_Dummy_Derivatives

The dummy derivative method as a follow up for pantelides algorithm to resolve singular systems. Is something like this planned?

[ChrisRackauckas]

Yes, definitely planned.

[MarkusLohmayer]

The following code does not work:

using ModelingToolkit
using DifferentialEquations
@variables t
D = Differential(t)
ps = @parameters mass₊m = 1.0 spring₊c = 10.0
sts = @variables mass₊p(t) spring₊q(t)
eqs = [
    0 ~ (-mass₊p) / mass₊m + D(spring₊q),
    0 ~ (-spring₊q) / spring₊c - D(mass₊p)
]
sys = ODESystem(eqs, t, sts, ps, name = :osc)
u0 = [
    spring₊q => 0.0,
    mass₊p => 10.0,
]
prob = ODAEProblem(structural_simplify(sys), u0, (0, 80.0))
sol = solve(prob, Tsit5(), abstol = 1e-8, reltol = 1e-8)

Without any Error, the "solution" is stuck at the initial condition!

In contrast, when defining eqs as

eqs = [
    D(spring₊q) ~ mass₊p / mass₊m,
    D(mass₊p) ~ (-spring₊q) / spring₊c
]

the behavior is as expected.

Chris mentioned on the Slack channel that this will be fixed by state selection, therefore I am posting to this issue.

The above code defines variables which look like namespace variables because the formerly defined eqs are actually the output of structural_simplify.

[ChrisRackauckas]

Partial state selection works nowadays