I’m interested in possibly using the DiffEqOperators functionality in DifferentialEquations.jl

I usually solve systems of nonlinear equations that are of the form

\frac{\partial A}{\partial z} = i\beta(z) A + F(A,z)

where A and \beta are complex arrays (i.e. \beta is a linear operator), and F is some nonlinear function of A. The way I do this is to transform the above into

And I have to transform back at each step to calculate F, and also transform back from G to A when saving the solution. Note that in general \beta depends on z, and so the operator needs to be updated as the integration progresses.

I currently solve the above by using the iterator interface to DifferentialEquations.jl, and mutating u after each step to reset it to A, before optionally saving and continuing. I get good performance from the Tsit5 algorithm, and I’m very happy (compared to when I was hand writing my own ODE solvers in C++).

However, I noticed the page on DiffEqOperators and wondered if that interface is designed to solve exactly this kind of problem?

Yes, it’s exactly for this kind of problem. More specifically, we have ExponentialUtilities.jl for doing the Krylov approximation to the matrix exponential

and then utilize this in the semilinear exponential integrators.

Can I ask if these methods are efficient for diagonal operators (which is my case)? For example is the diagonal trait from LinearMaps.jl used? Regardless, I will give this a try.

Are there any examples anywhere of how to pull all of the pieces together?

If you make an array operator where it’s a Diagonal matrix then it will make use of the overloads. There is more we should be doing like this. We don’t have such a trait to generalize it to matrix-free operations yet or anything like that.

There’s some random tests around

But we really need to get a whole tutorial series out for it. There’s still some more on the API we want to settle first though.