How to connect two ODEs?

[lindnemi]

I am trying to connect two oscillators, but don't manage to do it. I tried to generalize from the example given in the docs but get an error that only semi-explicit mass matrices are supported. Could you give me a hint?

using ModelingToolkit, OrdinaryDiffEq, Plots


### Attempt without edge

@parameters t P D c(t)
@variables θ(t) ω(t)
@derivatives Dt'~t

eqs1 = [Dt(θ) ~ ω,
        Dt(ω) ~ P - D * ω + sin(θ - c) ]

kura1 = ODESystem(eqs1, name=:kura1)


eqs2 = [Dt(θ) ~ ω,
        Dt(ω) ~ P - D * ω + sin(c - θ) ]


kura2 = ODESystem(eqs2, name=:kura2)



connections = [kura2.c ~ kura1.θ, kura1.c ~ kura2.θ]
sys = ODESystem(connections, t, [], [], systems=[kura1, kura2])


u0 = [kura1.θ => 0.,
      kura1.ω => 1.,
      kura2.θ => .5,
      kura2.ω => 0.3332]
    
p = [kura1.P => 1.,
     kura2.P => -1.,
     kura1.D => .1,
     kura2.D => .1]

tspan = (0., 10.)

prob = ODEProblem(sys, u0, tspan, p)
sol = solve(prob, Tsit5())
plot(sol)

[ChrisRackauckas]

Right now we don't do algebraic simplification, so the issue is as simple as

connections = [kura2.c ~ kura1.θ, kura1.c ~ kura2.θ]

needs to become

connections = [0 ~ kura1.θ - kura2.c, 0 ~ kura2.θ - kura1.c]

I'll turn this into the issue for tracking that.

(P.S. because this has a singular mass matrix you'd have to switch from Tsit5() to Rosenbrock23())

[lindnemi]

Thank you! Now ODEProblem produced an index error since it was attempting to acces u0 at index 6. Setting

u0 = [kura1.θ => 0.,
      kura1.ω => 1.,
      kura2.θ => .5,
      kura2.ω => 0.3332, 
      kura1.c => kura2.θ, 
      kura2.c => kura1.θ]

resolved the issue. However i think it is somewhat counter-intuitive to set initial conditions for parameters in the same array as for the dependent variables. I was expecting that the value of the coupling terms c would be determined by the connection equation. Anyway it just brought me to the next error by ODEProblem: ArgumentError: invalid index: nothing of type Nothing.

[ChrisRackauckas]

However i think it is somewhat counter-intuitive to set initial conditions for parameters in the same array as for the dependent variables. I was expecting that the value of the coupling terms c would be determined by the connection equation.

The c aren't parameters, they're algebraic variables.

[ChrisRackauckas]

This should all be good? Let me know if you have any more questions.

[lindnemi]

Unfortunately all is not good yet, i just didn't find the time to continue on this attempt.

At the moment I have

using ModelingToolkit, OrdinaryDiffEq, Plots

@parameters t P D c(t)
@variables θ(t) ω(t)
@derivatives Dt'~t

eqs1 = [Dt(θ) ~ ω,
       Dt(ω) ~ P - D * ω + sin(θ - c) ]

kura1 = ODESystem(eqs1, name=:kura1)


eqs2 = [Dt(θ) ~ ω,
       Dt(ω) ~ P - D * ω + sin(c - θ) ]

kura2 = ODESystem(eqs2, name=:kura2)


connections = [0 ~ kura2.c - kura1.θ; 0 ~ kura1.c - kura2.θ]
sys = ODESystem(connections, t, [kura1.c, kura2.c], [], systems=[kura1, kura2])


u0 = [kura1.θ => 0.,
      kura1.ω => 1.,
      kura2.θ => .5,
      kura2.ω => 0.3332,
      kura1.c => .5, #kura2.θ,
      kura2.c => .0] #kura1.θ]

p = [kura1.P => 1.,
     kura2.P => -1.,
     kura1.D => .1,
     kura2.D => .1,
     kura1.c => 0.,
     kura2.c => 0.]

tspan = (0., 10.)

prob = ODEProblem(sys, u0, tspan, p)
sol = solve(prob, Rosenbrock23())
plot(sol)

The construction of the ODEProblem works, but the integration terminates because of instabilities. What's wrong here?

And more fundamentally: Does MTK offer a way of connecting the two systems, without constructing a mass matrix? For what i want to achieve (coupling two oscillators) a mass matrix is not necessary, as c could be substituted with θ from the other system.

I am working with @FHell on NetworkDynamics.jl and wanted to look at MTK as a possible alternative for coupling heterogeneous systems.

[ChrisRackauckas]

And more fundamentally: Does MTK offer a way of connecting the two systems, without constructing a mass matrix? For what i want to achieve (coupling two oscillators) a mass matrix is not necessary, as c could be substituted with θ from the other system.

Not yet. That's the big issue #340 that I plan to get to soon.

[ChrisRackauckas]

The construction of the ODEProblem works, but the integration terminates because of instabilities. What's wrong here?

The index of this DAE isn't 1. Yeah, I think it just needs an explicit ODE form.