ODE of Vector of Matrices

New to Julia
,
1

I’m trying to construct/solve an ODE of several different sized matrices, but solving the ODE fails. The MWE below should demonstrate the issue:

using DifferentialEquations
using LinearAlgebra

function flow!(du,u,p,t)
    du[1]=u[1]
    du[2]=u[2]
end

UV=[Matrix(1.0I,3,3),Matrix(1.0I,5,5)]
timespan=(0.0,1.0)
prob=ODEProblem(flow!,UV,tspan)
sol=solve(prob)

ERROR: MethodError: no method matching zero(::Type{Matrix{Float64}})
Closest candidates are:
  zero(::Union{Type{P}, P}) where P<:Dates.Period at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/Dates/src/periods.jl:53
  zero(::DynamicPolynomials.Polynomial{C, T}) where {C, T} at /home/t248b327/.julia/packages/DynamicPolynomials/zw67w/src/poly.jl:26
  zero(::StatsBase.Histogram{T, N, E}) where {T, N, E} at /home/t248b327/.julia/packages/StatsBase/PGTj8/src/hist.jl:538
  ...
Stacktrace:
  [1] zero(x::Vector{Matrix{Float64}})
    @ Base ./abstractarray.jl:1085
...

Its seems to be trying to find an additive identity for the type Matrix{Float64}, but it doesn’t exist. I can run zero just fine on the individual matrices of the UV vector, but there is no general additive identity for a Matrix{Float64}.

What do I need to change type-wise to make this sort of formulation work (or alternative formulations if I’m completely misunderstanding how ODEproblem should work).

2

The problem is that the ODE solver needs to have an object in which all elements are indexed linearly, i.e. a Vector of matrices is not possible immediately.
There is a good work workaround though, provided by the package RecursiveArrayTools
https://diffeq.sciml.ai/stable/features/diffeq_arrays/#ArrayPartitions

using RecursiveArrayTools
A = ArrayPartition(Matrix(1.0I,3,3),Matrix(1.0I,5,5))

where A.x[1] is your first matrix and so on.
this partition array is just a big vector under the hood and therefore you can use it for ODE’s.

1 Like
3

Okay, that makes more sense. So the ArrayPartition basically just puts it in a form where the ODE solver will accept it, but also gives me special accessors so I can still get to my original matrices.

It seems to work just fine using the ArrayPartition. Thank you for the help!

4

Exactly. Maybe I should be a bit more precise about why this is necessary: Any ODE solver at has to somehow add together the derivative with the state object, in a simple Euler scheme we could implement ourselves there might be something like

for i in eachindex(Deriv)
  State[i] += Deriv[i] * deltaT
end

Essentially, the ODE solver assumes that there are a few functions being defined for the types of your vector, such as +/-, norm (to compute errors in adaptive solvers) and apparently also zero, I assume to initialize one or several derivative vectors with zeros of the type you use. However, in your previous case we had for instance that A[1] is a matrix. For a matrix, there was no defined function “zero”, so an error was thrown.
In principle, there are two workarounds for this: You could implement all the functions that are needed for your type (Matrix), or we throw all elements together in a big array, such that its elements are numbers, for which zero() and the other functions are already defined.
RecursiveArrayTools does the latter, so that you can get all elements of both your matrices by iterating over A.
The clever bit is just that it also allows you to recover your matrices for convenience.
Hope that clears it up a bit, I remember being really impressed by this module when I found it, after trying my own horrible ideas :slight_smile: