Defing ODEsystem in function and initial conditions outside the function in ModelingToolkit.jl

Specific Domains Numerics
1

Hi,

I try to initialize an ODE system in a function and then assign initial conditions and parameters through any kind of named structure (dict, vector(tuples), or whatever). However, the only way to set parameters outside the function is then through arrays, which is a source of error when they get large…
I’m new to Julia and try to find a good way to solve this that scales well for ensembles of ODEs later.
I tried this:

using ModelingToolkit, DifferentialEquations
function init_ODE()
        @parameters g
        @variables t, x(t), y(t), c_gp(t)

        D = Differential(t)
        particle_equations = [
                c_gp    ~ - x./ g,
                D(x)   ~  c_gp ,
                D(y)   ~  x-  u(t)
                ]
        return particle_equations
end

@register_symbolic u(t)
@named fol_separate = ODESystem(init_ODE())
fs = structural_simplify(fol_separate)

# external forcing
u(t) = 7 + 2.0 .* sin.(t .* π ./(60*30))
z0 = [ x => 0.0,  y => 0.0]
params = [ g  => 9 ]

prob = ODEProblem(fs, z0, (0.0,  60*60*5) , params) #(0.0,60.0*60.0*10)
sol  = solve(prob, saveat=60*5)

gives an UndevVarErr for x, y, and g, which makes sense, they were never defined outside the function.
If, alternatively, I do:

z0 = [ fs.x => 0.0,  fs.y => 0.0]
params = [ fs.g  => 9 ]

The ODEProblem returns:

ArgumentError: Term{Real, Base.ImmutableDict{DataType, Any}}[x(t), y(t)] are missing from the variable map.

The only way is to write z0, params as arrays since @variables are only defined in the function.
Does someone have an idea of how to initialize a generic ODE problem with the ModelingToolkit properly and then define the parameters in a named structure?

In the end, I want a large ensemble of ODE systems that only differ by parameters and forcing. containing the ODEs and algebraic expressions in a function would be very useful.

2

You can

@unpack x,y = fs
z0 = [ x => 0.0,  y => 0.0]
...

to load them into the current scope. Or just declare them again in the scope you want to create the initial condition using @variables again

@variables t, x(t), y(t)
z0 = [ x => 0.0,  y => 0.0]
3

Thanks! that does it. repeating the declaration is another source of error, but even when I @unpack I have to know what names to unpack, but thats okey.
thanks!
m

4

I agree it would be nice if you could just use fs.x instead. I think there are open issues about this on ModelingToolkit’s repo, or at least there have been in the past, but perhaps there are reasons that allowing it could be problematic.

5

Any further development on this @mochell ? I’m currently running into the same issue in a similar use case. If not, I will try redeclaring the variables outside of the function scope.

6

Hey,

I should have followed up. I figured out a workaround by declaring the same variables in the function and the executing script. Likely not ideal but works fine for now:

using ModelingToolkit, DifferentialEquations

function init_vars()
        @parameters g
        @variables t, x(t), y(t), c_gp(t)

        @parameters r_g C_α C_φ g C_e
        return  t, x, y, c_gp, g
end 

function init_ODE(u)
        t, x, y, c_gp, g = init_vars()

        D = Differential(t)
        particle_equations = [
                c_gp    ~ - x./ g,
                D(x)   ~  c_gp ,
                D(y)   ~  x-  u(t)
                ]
        return particle_equations
end


u(t) = 7 + 2.0 .* sin.(t .* π ./(60*30))

t, x, y, c_gp, g = init_vars()
z0 = [ x => 0.0,  y => 0.0]
params = [ g  => 9 ]

@named particle_system0 = ODESystem( init_ODE(u) )
particle_system1 = structural_simplify(particle_system0)
problem    = ODEProblem(particle_system1, z0, (0.0,  60*60) , params)

sol  = solve(problem, saveat=60*5)

Both functions can then be put in a module and called from another script.