How to tell ModelingToolkit that two variables have the same initial conditions without computing them

Specific Domains Modelling & Simulations
,
1

Suppose that we have a differential-algebraic system of equations, with two variables f(t) and g(t) such that f(0) = g(0). Suppose moreover that g(t) is algebraic and g(0) could be explicitly calculated.
Is it possible to avoid the latter calculation and let MTK figure it out by itself?

Here is a MWE of what we are trying to do, which returns an error:

using ModelingToolkit, DifferentialEquations

@variables t
D = Differential(t)

function system_f(; name)
    @variables f(t) g(t)
    eqs = [D(f) ~ g - f]
    ODESystem(eqs; name)
end

function system_g(; name)
    @variables g(t)
    eqs = [g ~ t]
    ODESystem(eqs; name)
end

@named f = system_f()
@named g = system_g()

connected_eqs = [
    f.g ~ g.g,
]

@named _model = ODESystem(connected_eqs, t)
@named model = compose(_model, [f, g])

model = structural_simplify(model)
prob = ODEProblem(model, [f.f => g.g], (0.0, 10.0))
sol = solve(prob)

In the above example, we have that g(t)=t, thus f(0)=g(0)=0.
When g(t) is much more complicated, is there a way to avoid having to provide g(0) to MTK explicitly?

2 Likes
2

Define a parameter and make both initial conditions equal to that parameter.

3

Hi Chris, thanks for your answer.

Regarding your proposal, how can we set the parameter value to be dependent on g(t_0)?

The initial value of g is automatically computed and dependent on t_0 (from the timespan passed to ODEProblem), thus we cannot assign a constant value to it.
Since we don’t need to provide an initial value of g, we also don’t want to provide an initial value to f (or to a parameter p), because it can be clearly derived from the problem.

4

One -quite inelegant- way I’m managing to address the problem is to solve the system 2 times, where the first time is merely used to get the initial value of g, as follows:

  • first set both initial conditions equal to an arbitrary value
  • solve the system, getting the solution pre_sol,
  • solve again the system, this time with initial condition f(0)=\text{pre_sol}.g(0)

For example, here are the lines to substitute to the last two ones in the initial MWE:

pre_prob = ODEProblem(model, [f.f=>42,g.g=>42], (0.0, 10.0)) 
pre_sol = solve(pre_prob)
prob = ODEProblem(model, [f.f=>pre_sol[g.g, 1],g.g=>pre_sol[g.g, 1]], (0.0, 10.0))
sol = solve(prob)

Of course, I suppose that there’s a more direct solution.

5

Open an issue.

6

Opened: Solving a system with initial condition of the form $f(0) = g(0)$ · Issue #1911 · SciML/ModelingToolkit.jl · GitHub

7

Just a small update: we added a simpler MWE in the MTK issue.
We’ve been hitting this limitation again lately, in trying to port some famous models from Vensim to a Julia library.

3 Likes