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[1] = u[1] - u[2][1]
    du[2] .= -5 * u[2]
    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!