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?
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.