# How can I define a vector in ModelingToolkit

**URL:** <https://discourse.julialang.org/t/how-can-i-define-a-vector-in-modelingtoolkit/107278>\
**Category:** Modelling & Simulations\
**Tags:** question, package, modelingtoolkit\
**Created:** [December 7, 2023, 3:13pm UTC](https://discourse.julialang.org/t/how-can-i-define-a-vector-in-modelingtoolkit/107278 "2023-12-07T15:13:00Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [December 7, 2023, 3:13pm UTC](https://discourse.julialang.org/t/how-can-i-define-a-vector-in-modelingtoolkit/107278/1 "2023-12-07T15:13:00Z")

</div>

The following code fails:

```julia
# Example one:
# Falling mass.

using ModelingToolkit, OrdinaryDiffEq

G_EARTH = Float64[0.0, 0.0, -9.81] # gravitational acceleration

# model
@variables t pos(t)[1:3]=[0.0, 0.0, 0.0]
@variables vel(t)[1:3]=[0.0, 0.0, 50.0] 
@variables acc(t)[1:3]=[0.0, 0.0, -9.81] 
D = Differential(t)

eqs = [D(pos) ~ vel,
       D(vel) ~ acc,
       acc ~ G_EARTH]

@named sys = ODESystem(eqs, t)

```

with the message:

```julia
ERROR: LoadError: axes of Differential(t)((pos(t))[1:3]) not known
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] axes(A::Symbolics.Arr{Num, 1})
   @ Symbolics ~/.julia/packages/Symbolics/5BibL/src/arrays.jl:537 [inlined]
 [3] combine_axes
   @ ./broadcast.jl:524 [inlined]
 [4] instantiate
   @ ./broadcast.jl:306 [inlined]
 [5] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}, Nothing, typeof(~), Tuple{Num, Num}})
   @ Base.Broadcast ./broadcast.jl:903 [inlined]
 [6] ~(lhs::SymbolicUtils.BasicSymbolic{Symbolics.Arr{Num, 1}}, rhs::Symbolics.Arr{Num, 1})
   @ Symbolics ~/.julia/packages/Symbolics/5BibL/src/equations.jl:104
 [7] top-level scope
   @ ~/repos/Tether/Tether_01.jl:17
 [8] include(fname::String)
   @ Base.MainInclude ./client.jl:489
 [9] top-level scope
   @ REPL[1]:1
in expression starting at /home/ufechner/repos/Tether/Tether_01.jl:17

julia> include("Tether_01.jl")
ERROR: LoadError: axes of Differential(t)((pos(t))[1:3]) not known
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] axes(A::Symbolics.Arr{Num, 1})
   @ Symbolics ~/.julia/packages/Symbolics/5BibL/src/arrays.jl:537 [inlined]
 [3] combine_axes
   @ ./broadcast.jl:524 [inlined]
 [4] instantiate
   @ ./broadcast.jl:306 [inlined]
 [5] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}, Nothing, typeof(~), Tuple{Num, Num}})
   @ Base.Broadcast ./broadcast.jl:903 [inlined]
 [6] ~(lhs::SymbolicUtils.BasicSymbolic{Symbolics.Arr{Num, 1}}, rhs::Symbolics.Arr{Num, 1})
   @ Symbolics ~/.julia/packages/Symbolics/5BibL/src/equations.jl:104
 [7] top-level scope
   @ ~/repos/Tether/Tether_01.jl:14
 [8] include(fname::String)
   @ Base.MainInclude ./client.jl:489
 [9] top-level scope
   @ REPL[1]:1
in expression starting at /home/ufechner/repos/Tether/Tether_01.jl:14

```

What am I doing wrong?

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [December 7, 2023, 3:41pm UTC](https://discourse.julialang.org/t/how-can-i-define-a-vector-in-modelingtoolkit/107278/2 "2023-12-07T15:41:30Z")

</div>

This does not give a syntax error:

```julia
# Example one:
# Falling mass.

using ModelingToolkit, OrdinaryDiffEq

G_EARTH = Float64[0.0, 0.0, -9.81] # gravitational acceleration

# model
@variables t pos(t)[1:3]=[0.0, 0.0, 0.0]
@variables vel(t)[1:3]=[0.0, 0.0, 50.0] 
@variables acc(t)[1:3]=[0.0, 0.0, -9.81] 
D = Differential(t)

eqs = vcat(D.(pos) ~ vel,
           D.(vel) ~ acc,
           acc .~ G_EARTH)

@named sys = ODESystem(eqs, t)

```

Lets see if it works…

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [December 7, 2023, 5:35pm UTC](https://discourse.julialang.org/t/how-can-i-define-a-vector-in-modelingtoolkit/107278/3 "2023-12-07T17:35:43Z")

</div>

Works fine:

```julia
# Example one: Falling mass.
using ModelingToolkit, OrdinaryDiffEq, PyPlot

G_EARTH = Float64[0.0, 0.0, -9.81] # gravitational acceleration

# model
@variables t pos(t)[1:3]=[0.0, 0.0, 0.0]
@variables vel(t)[1:3]=[0.0, 0.0, 50.0] 
@variables acc(t)[1:3]=[0.0, 0.0, -9.81] 
D = Differential(t)

eqs = vcat(D.(pos) ~ vel,
           D.(vel) ~ acc,
           acc .~ G_EARTH)

@named sys = ODESystem(eqs, t)
simple_sys = structural_simplify(sys)

duration = 10.0
dt = 0.02
tspan = (0.0, duration)
ts = 0:dt:duration
u0 = zeros(6)
u0[6] = 50.0

prob = ODEProblem(simple_sys, u0, tspan)
sol = solve(prob, Rodas5(), dt=dt, saveat=ts)

X = sol.t
POS_Z = sol(X, idxs=pos[3])
VEL_Z = sol(X, idxs=vel[3])

plot(X, POS_Z, color="green")
xlabel("time [s]")
grid(true)
twinx()
ylabel("vel_z [m/s]") 
plot(X, VEL_Z, color="red") 
nothing

```

Only remaining question:  
Why does `POS = sol(X, idxs=pos[1:3])` not work?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [December 7, 2023, 6:27pm UTC](https://discourse.julialang.org/t/how-can-i-define-a-vector-in-modelingtoolkit/107278/4 "2023-12-07T18:27:13Z")

</div>

> [@ufechner7](#):
>
> Why does `POS = sol(X, idxs=pos[1:3])` not work?

Symbolic arrays have a number of problems, you might need to do `idxs = collect(pos)` for now

---

<div class="post-metadata">

**Author:** ![jondavis847](https://avatars.discourse-cdn.com/v4/letter/j/da6949/32.png) [@jondavis847](https://discourse.julialang.org/u/jondavis847)\
**Post date:** [December 7, 2023, 8:10pm UTC](https://discourse.julialang.org/t/how-can-i-define-a-vector-in-modelingtoolkit/107278/5 "2023-12-07T20:10:10Z")

</div>

you can do

```julia

julia> sol[pos]
501-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.0]
 [0.0, 0.0, 0.9980380000000015]
 [0.0, 0.0, 1.9921519999999922]
 [0.0, 0.0, 2.982341999999982]
 [0.0, 0.0, 3.96860799999997]
 [0.0, 0.0, 4.950949999999957]
 ⋮
 [0.0, 0.0, 14.260949999999628]
 [0.0, 0.0, 13.316607999999672]
 [0.0, 0.0, 12.368341999999714]
 [0.0, 0.0, 11.416151999999672]
 [0.0, 0.0, 10.460037999999674]
 [0.0, 0.0, 9.499999999999718]

```

---

<div class="post-metadata">

**Author:** ![RobbesU](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/robbesu/32/212459_2.png) [@RobbesU](https://discourse.julialang.org/u/RobbesU)\
**Post date:** [January 29, 2025, 11:12am UTC](https://discourse.julialang.org/t/how-can-i-define-a-vector-in-modelingtoolkit/107278/6 "2025-01-29T11:12:55Z")

</div>

Above code block doesn’t work anymore under `ModelingToolkit v9.61.0` (Julia 1.11.3) but works with a few adjustments.

```julia
duration = 5.0
tspan = (0.0, duration)
dt = 0.02
ts = 0:dt:duration
	
G_EARTH = Float64[0.0, -0.0, -9.81] # gravitational acceleration
G_MOON = G_EARTH * 0.166
	
# model
@variables (pos(t))[1:3] = [0.0, 0.0, 0.0] # defaults to 'at rest at origin'
@variables (vel(t))[1:3] = [0.0, 0.0, 0.0] 
@variables (acc(t))[1:3] # NO defaults as balanced by force!

@parameters (force)[1:3] # do not provide default or structural_simplify will remove eqs where f=0 !
	
eqs = [
	collect( D.(pos) .~ vel )... # .~ instead of just ~ !
	collect( D.(vel) .~ acc )...
    collect( acc .~ force )...
]
	
@named sys = ODESystem(eqs, t)
simple_sys = structural_simplify(sys)
	
u0 = [pos=>[0.,0.,0.], vel=>[0.,0.,+10.0] ] # MTK unknows are pos & vel: unknowns(simple_sys) 
p0 = [[0.,0.,-9.81]] # [G_EARTH] 

prob = ODEProblem(simple_sys, u0, tspan, p0)
sol = solve(prob, Rodas5(), dt=dt, saveat=ts)
plot(sol)

```

When constructing as a component block it becomes as follows (I added a force as well as default gravity). Note the requirement to state the initial condition keys as symbolic `:pos` and `:vel` if using the Dict version of `u0`:

```julia
# model as component
@component function TestMass3D(;name, m, g=G_EARTH)
	
	vars = @variables begin
		(pos(t))[1:3] = [0.0, 0.0, 0.0] # defaults to 'at rest at origin'
		(vel(t))[1:3] = [0.0, 0.0, 0.0] 
		(acc(t))[1:3] # NO defaults as balanced by force!
	end
		
	pars = @parameters (force)[1:3] # do not provide default or structural_simplify will remove eqs where f=0 !
	
	eqs = [
		collect( D.(pos) .~ vel )...
		collect( D.(vel) .~ acc )...
	    collect( acc .~ force./m .+ g )...
	]

	return ODESystem(eqs, t, vars, pars;name)
end

@named sys = TestMass3D(m=1, g= G_MOON)
simple_sys = structural_simplify(sys)
	
u0 = [:pos=>[0.,0.,0.], :vel=>[0.,0.,+10.0] ] # MTK unknows are pos & vel: check with unknowns(simple_sys) 
p0 = [[0.,0.,+2.] ] # note vector inside array as only 1 parameter 

prob = ODEProblem(simple_sys, u0, tspan, p0)

sol = solve(prob, Rodas5(), dt=dt, saveat=ts)
plot(sol)

```

(this also links with [https://discourse.julialang.org/t/broadcast-inside-of-mtkmodel-macro-causes-precompilation-failures-when-put-into-a-package/107277/18](https://discourse.julialang.org/t/broadcast-inside-of-mtkmodel-macro-causes-precompilation-failures-when-put-into-a-package/107277/18)