Best way to solve system of ODEs with different types in DifferentialEquations.jl

I have two fields, `A` (complex scalar field) and `B` (real vector field), which are sampled on staggered grids. This means that (in, say, three dimensions) `A` is an complex array with N^3 elements, while `B` is a real array with 3*(N-1)^3 elements. I also have corresponding equations of motion with `dA = fA(A,B)` and `dB = fB(A,B)`. I’m rolling my own integrator, but I would prefer to switch to DifferentialEquations.jl.

I can see a trivial, even if inconvenient, solution to the difference of array shapes by flattening and stacking both into a vector and then `view`-ing `reshape`-ing them into correct form inside the problem definition.

My bigger concern is the complex/real split between the arrays. This type difference has physical relevance and is enforced through the code base, so I would prefer to not cast the real field into a complex-valued field just for sake of this numerical step. There is also space, where in practical situations, N will be large (easily N=10^4, could increase magnitude or two) and I’m not too keen on almost doubling the memory footprint just to store billions of zeros.

What is the most graceful way of solving a system of equations like this?

1 Like

An ArrayPartition should be a good data structure for this.

2 Likes