Solving a system of matrix differential equations using DifferentialEquations.jl

General Usage
,
1

How shall I initialize my initial condition for the system of matrix differential equations using DifferentialEquations.jl? As a simple example, suppose I want to solve the following trivial system of equations:

function g!(du,u,p,t)
   du[1] .= u[1] + ones(3,3)
   du[2] .= ones(3,3)
end

If I then define my problem as

u0 = [zeros(3,3),zeros(3,3)]
prob = ODEProblem(g!,u0,(0,1.0))

and try to solve it solve(prob,Tsit5()), I get the error:

ERROR: MethodError: no method matching zero(::Type{Matrix{Float64}})
Closest candidates are:
zero(::Union{Type{P}, P}) where P<:Dates.Period at C:\ …\Julia-1.8.5\share\julia\stdlib\v1.8\Dates\src\periods.jl:53
zero(::AbstractIrrational) at irrationals.jl:150
zero(::FillArrays.AbstractOnes{T, N}) where {T, N} at C:\ ….julia\packages\FillArrays\sIOcj\src\FillArrays.jl:605
…

I also tried to pass my initial condition as a tuple of matrices as well as as a static array with an out-of-place definition of the problem, but nothing worked. Are systems of matrix equations supported by the package, and if so, how do I correctly define the problem?

2

Your state array has to be an abstractArray containing numbers, yours is an array of arrays.

You can either include your two matrices into a higherdimensional array u0 = zeros(2,3,3) , or you use more specialized types such as from RecursiveArrayTools.jl

3

How could I not think of that… So this is what I ended up with

function g!(du,u,p,t)
   du[1,:,:] = u[1,:,:] + ones(3,3)
   du[2,:,:] = ones(3,3)
end

which seems to do the trick.

4

See the new error message:

ERROR: Non-Number element type inside of an Array detected.
Arrays with non-number element types, such as
Array{Array{Float64}}, are not supported by the
solvers.

If you are trying to use an array of arrays structure,
look at the tools in RecursiveArrayTools.jl. For example:

If this was a mistake, promote the element types to be
all the same. If this was intentional, for example,
using Unitful.jl with different unit values, then use
an array type which has fast broadcast support for
heterogeneous values such as the ArrayPartition
from RecursiveArrayTools.jl. For example:

using RecursiveArrayTools
u0 = ArrayPartition([1.0,2.0],[3.0,4.0])
u0 = VectorOfArray([1.0,2.0],[3.0,4.0])

are both initial conditions which would be compatible with
the solvers.

Element type:
Matrix{Float64}

5

Thanks! Is using RecusiveArrayTools.jl preferable compared to using higher-dimensional arrays?

6

no. using a higher dimensions array is preferred for fixed sizes.

1 Like