Matrix coupled equations

[jebej]

I have a system of many coupled equations, where the output of each equation is a matrix. When you have a coupled system of scalar functions, you can just set up u and du to be Vectors, and unpack the system in a single vector ODE. What can you do when each equation instead outputs a matrix?

I tried simply choosing my u to be a vector of matrices, but the solver fails while trying to set the cache up:

MethodError: no method matching zero(::Type{Array{Complex{Float64},2}})
in #solve#1419 at OrdinaryDiffEq\src\solve.jl:6
in  at base\<missing>
in #init#1420 at OrdinaryDiffEq\src\solve.jl:246
in alg_cache at OrdinaryDiffEq\src\caches\verner_caches.jl:95

I also tried putting everything in a 3D array, and using views, but that didn't work either.

[jebej]

hmm things work if you concatenate the matrices horizontally, which I assume is going to be the most efficient once views stop allocating.

[ChrisRackauckas]

I don't quite understand your issue. If you use u0 as a matrix/3-tensor though then it will use that in the system. Can I get a minimal example of the problem?

[jebej]

Thanks Chris, I can't reproduce the problem, so I must have made a mistake (probably when trying out different function styles).

Nevertheless, concatenating the matrices (whether in the 2nd or 3rd dimension) is not ideal due to all the views that need to be made for each function execution (needed for matrix multiplication). Could the solver somehow be made to work with a vector of matrices instead?

[ChrisRackauckas]

Could the solver somehow be made to work with a vector of matrices instead?

It would take a lot of work to make arrays of arrays work. That's #254 . Arrays of static arrays work, so a vector of static matrices will do the trick if those matrices are small enough. Another way to do this is using an ArrayPartition from RecursiveArrayTools.jl. This will only work with Julia-defined solvers. This will do well on Julia v0.7, but on Julia v0.6 there's an issue with broadcast and so it will be slow on these methods: SciML/OrdinaryDiffEq.jl#106. When Julia v0.7 comes out though, it will be a fantastic option.

Here's a blog post demonstrating array partitions: http://www.stochasticlifestyle.com/solving-systems-stochastic-pdes-using-gpus-julia/

[jebej]

RecursiveArrayTools seems to use a tuple of arrays as storage, is it possible to use that directly? Maybe by manually specifying some needed functions?

[ChrisRackauckas]

You'd have to define a recursive broadcast, which is the most difficult part of the code.

[jebej]

Understood, thanks Chris.