Possible Inconsistent Initial Values of Some State Variables

[ggarza]

In the following, the values of the state variables ra=r_(t) and sa=s_(t) can be determined for any initial value, yet MTK still identifies them as state variables that must be given an initial value. And giving them an initial value other than r_(0) and s_(0) leads to issues with the solution.

using ModelingToolkit
using DifferentialEquations

## Fox-Rabbit Pursuit
@variables  t x(t) y(t) R(t) ra(t) sa(t) dra(t) dsa(t)
@parameters k
D = Differential(t)

# rabbit's path, an outward spiral
r_(t) = sqrt(1+t) * cos(t)
s_(t) = sqrt(1+t) * sin(t)

@named rabbitfox_model = ODESystem([
    ra ~ r_(t),
    sa ~ s_(t),
    D(ra) ~ dra,
    D(sa) ~ dsa,
    D(x) ~ R * (ra - x),
    D(y) ~ R * (sa - y),
    R ~ k * sqrt(dra^2 + dsa^2)/
        sqrt((ra - x)^2 + (sa - y)^2),
])

prob = ODEProblem(structural_simplify(rabbitfox_model),
            [x=>3, y=>0, ra=>r_(0), sa=>s_(0)],
            (0, 5.071),
            [k=>1.1]
            )

sol = solve(prob)

(As a newbie to MTK), I originally thought that the initial values for ra and sa would be automagically determined.

This issue is in response to this discussion: https://discourse.julialang.org/t/argumenterror-in-fox-rabbit-pursuit-problem-with-modelingtoolkit-jl/61222/2

[ChrisRackauckas]

I don't understand the question: this works for me?

[ggarza]

Yes it works.

My issue is that given the original problem it only makes sense for the initial values of ra and sa to be r_(0) and s_(0) but MTK is allowing us to set other initial values. It seems to me that MTK should detect that the only initial values actually needed are for x and y.

That is, having

prob = ODEProblem(structural_simplify(rabbitfox_model),
            [x=>3, y=>0],
            (0, 5.071),
            [k=>1.1]
            )

should be sufficient to solve the problem.

I'm not sure if this is the bug that @YingboMa thought he saw. Honestly, I don't really know how MTK works yet.

[ChrisRackauckas]

Oh I see, because ra and sa are algebraic. But then D(ra) ~ dra and D(sa) ~ dsa don't make any sense: it's already known, and the derivative isn't constant, so delete that.

using ModelingToolkit
using DifferentialEquations

## Fox-Rabbit Pursuit
@variables  t x(t) y(t) R(t) ra(t) sa(t) dra(t) dsa(t)
@parameters k
D = Differential(t)

# rabbit's path, an outward spiral
r_(t) = sqrt(1+t) * cos(t)
s_(t) = sqrt(1+t) * sin(t)

using ForwardDiff
dr_(t) = ForwardDiff.derivative(r_,t)
ds_(t) = ForwardDiff.derivative(s_,t)


@named rabbitfox_model = ODESystem([
    ra ~ r_(t),
    sa ~ s_(t),
    D(x) ~ R * (ra - x),
    D(y) ~ R * (sa - y),
    R ~ k * sqrt(dr_(t)^2 + ds_(t)^2)/
        sqrt((ra - x)^2 + (sa - y)^2),
])

prob = ODEProblem(structural_simplify(rabbitfox_model),
            [x=>3, y=>0],
            (0, 5.071),
            [k=>1.1]
            )

sol = solve(prob)

seems to work fine. I'd put making that more automatic as a low priority issue since it's a pretty special case with an easy user-side fix.

[ggarza]

ok. Thank you. Feel free to close this issue.