# ODE with heterogeneous input through RecursiveArrayTools

Hi all
I am trying to solve a first-order ODE with initial conditions input as a heterogeneous array using of RecursiveArrayTools.jl https://docs.sciml.ai/DiffEqDocs/stable/features/diffeq_arrays/#ArrayPartitions.

I cannot figure out what causes this.

using OrdinaryDiffEq, LinearAlgebra, ForwardDiff, RecursiveArrayTools

A    = 0.225
σ_oe = 490e6
C    = 1.25
K    = 2.65e-5
L    = 14.4

# Inline functions
s(x) = x - 1.0/3.0*tr(x)*diagm(ones(size(x,1)))
σ_eff(x1, x2) = sqrt(3/2)*norm(s(x1)-x2,2)
β(x1, x2) = 1/σ_oe * (σ_eff(x1, x2) + A*tr(x1)-σ_oe)

function CD_ODE!(du,u,p,t)
α .= u.x[1]
D .= u.x[2]
σ .= p[1](t)
dσ .= p[2](t)

dβ = 1/(σ_oe + C*σ_eff(σ, α))*(3/2*tr((s(σ).-α)*dσ')/σ_eff(σ, α) + A*tr(dσ))

if  dβ > 0.0 && β(σ, α) > 0.0
du.x[1] .= dβ*C*(s(σ).-α)
du.x[2] .= dβ*K*exp(L*β(σ, α))
else
du.x[1] .= 0.0
du.x[2] .= 0.0
end
end

α₀ = zeros(Float64, 3,3)
D₀ = 0.0::Float64
#u₀ = [α₀ D₀]
u₀ = ArrayPartition(α₀, D₀)

t_end = 100.
tspan = (0.0,t_end)

f = 1.0

signal  = t -> 500e6*cos(2*pi*f*t)*diagm(ones(size(α₀,1)))
dsignal = t -> ForwardDiff.derivative(signal, t)

# p = [signal]
p = [signal dsignal]

prob = ODEProblem(CD_ODE!,u₀,tspan,p)
sol = solve(prob, Tsit5())

which gives me the following error:

ERROR: MethodError: no method matching similar(::Float64, ::Type{Float64})
Closest candidates are:
similar(::Union{Adjoint{T, var"#s886"}, Transpose{T, var"#s886"}} where {T, var"#s886"<:(AbstractVector)}, ::Type{T}) where T at ~/.julia/juliaup/julia-1.8.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.8/LinearAlgebra/src/adjtrans.jl:207
similar(::Union{Adjoint{T, S}, Transpose{T, S}} where {T, S}, ::Type{T}) where T at ~/.julia/juliaup/julia-1.8.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.8/LinearAlgebra/src/adjtrans.jl:211
similar(::Union{Adjoint{T, S}, Transpose{T, S}} where {T, S}, ::Type{T}, ::Tuple{Vararg{Int64, N}}) where {T, N} at ~/.julia/juliaup/julia-1.8.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.8/LinearAlgebra/src/adjtrans.jl:212
...
Stacktrace:
[1] (::RecursiveArrayTools.var"#52#53"{Float64, ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}})(i::Int64)
@ RecursiveArrayTools ~/.julia/packages/RecursiveArrayTools/9mtCv/src/array_partition.jl:59
[2] ntuple
@ ./ntuple.jl:49 [inlined]
[3] ArrayPartition
@ ~/.julia/packages/RecursiveArrayTools/9mtCv/src/array_partition.jl:31 [inlined]
[4] similar
@ ~/.julia/packages/RecursiveArrayTools/9mtCv/src/array_partition.jl:59 [inlined]
[5] alg_cache(alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, u::ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, rate_prototype::ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, #unused#::Type{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, uprev::ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, uprev2::ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, f::ODEFunction{true, SciMLBase.AutoSpecialize, typeof(CD_ODE!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, t::Float64, dt::Float64, reltol::Float64, p::Matrix{Function}, calck::Bool, #unused#::Val{true})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vIqLB/src/caches/low_order_rk_caches.jl:522
[6] __init(prob::ODEProblem{ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, Tuple{Float64, Float64}, true, Matrix{Function}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(CD_ODE!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vIqLB/src/solve.jl:322
[7] __init (repeats 5 times)
@ ~/.julia/packages/OrdinaryDiffEq/vIqLB/src/solve.jl:10 [inlined]
[8] __solve(::ODEProblem{ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, Tuple{Float64, Float64}, true, Matrix{Function}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(CD_ODE!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vIqLB/src/solve.jl:5
[9] __solve
@ ~/.julia/packages/OrdinaryDiffEq/vIqLB/src/solve.jl:1 [inlined]
[10] #solve_call#22
@ ~/.julia/packages/DiffEqBase/jb48J/src/solve.jl:496 [inlined]
[11] solve_call(_prob::ODEProblem{ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, Tuple{Float64, Float64}, true, Matrix{Function}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(CD_ODE!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False})
@ DiffEqBase ~/.julia/packages/DiffEqBase/jb48J/src/solve.jl:466
[12] solve_up(prob::ODEProblem{ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, Tuple{Float64, Float64}, true, Matrix{Function}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(CD_ODE!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, p::Matrix{Function}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ DiffEqBase ~/.julia/packages/DiffEqBase/jb48J/src/solve.jl:917
[13] solve_up
@ ~/.julia/packages/DiffEqBase/jb48J/src/solve.jl:890 [inlined]
[14] solve(prob::ODEProblem{ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, Tuple{Float64, Float64}, true, Matrix{Function}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(CD_ODE!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ DiffEqBase ~/.julia/packages/DiffEqBase/jb48J/src/solve.jl:827
[15] solve(prob::ODEProblem{ArrayPartition{Float64, Tuple{Matrix{Float64}, Float64}}, Tuple{Float64, Float64}, true, Matrix{Function}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(CD_ODE!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False})
@ DiffEqBase ~/.julia/packages/DiffEqBase/jb48J/src/solve.jl:817
[16] top-level scope
@ ~/UserDir/MWE.jl:50

Any ideas for what’s going wrong here? Thanks!

Try D₀ = [0.0]? Having one part just be a scalar might be the issue. I think they all need to be arrays if you want to mutate it.

1 Like

Thanks @ChrisRackauckas! Now it solves (…but doesn’t enter the non-trivial part of the conditions, but that’s for me to debug ).

I also modified these lines within the ODE function in the MWE:

α = u.x[1]
D = u.x[2]
σ = p[1](t)
dσ = p[2](t)