 # Nested vector variables in ODE problem definitions

I am trying to define an ODE in which the state variables are nested in arrays. Here is a toy example:

``````using DifferentialEquations

function f!(du, u, p, t)
du = u - u
du .= -5 * u
return nothing
end

u₀ = [1.0, [2.0, 3.0]]
tspan = (0.0, 10.0)

prob = ODEProblem(f!, u₀, tspan)

sol = solve(prob, Tsit5())
``````

I assume this type of thing is supported since the `solve` function makes a call to the RecursiveArrayTools function `recursive_unitless_bottom_eltype` to get the base element type of the nested array. But here it throws a stack overflow error because the methods of `recursive_unitless_bottom_eltype` look like:

``````recursive_unitless_bottom_eltype(a) = recursive_unitless_bottom_eltype(typeof(a))
recursive_unitless_bottom_eltype(a::Type{T}) where T = recursive_unitless_bottom_eltype(eltype(a))
recursive_unitless_bottom_eltype(a::Type{T}) where {T<:AbstractArray} = recursive_unitless_bottom_eltype(eltype(a))
recursive_unitless_bottom_eltype(a::Type{T}) where {T<:Number} = eltype(a) == Number ? Float64 : typeof(one(eltype(a)))
``````

Since the `typeof(u₀)` is `Array{Any,1}`, the second pass goes into the second method. Since `eltype(Array{Any,1})` is `Any` and `eltype(Any)` is `Any`, we recurse indefinitely.

Is my use case supported or am I doing something that just shouldn’t be done?

This is one thing I haven’t been able to solve without a type because you cannot declare broadcast as recursive. `u₀ = VectorOfArray([[1.0], [2.0, 3.0]])` works though, and so does `ArrayPartition`.

2 Likes

Oh, cool, I’ll have to give that a try. Thanks!