Hey ya’ll, for the coupled channels scattering RBM stuff we’re having some difficulties using scipy.integrate for a coupled second order linear system, where each sub-matrix in the block matrix is sparse. (This corresponds to the Schroedinger eqn w/ coupling between channels).
I think the issue is that scipy is set up for first order systems, so for an ND 2nd order problem, you have to convert it to a 2ND 1st order problem. Beyond that, it just uses Runge-Kutta and iterates, which doesn’t take advantage of that fact that the system is linear and the sub matrices are sparse.
The end result is we’re getting some weird numerical instability and oscillations in the resulting wavefunctions, This could definitely be user error, but if anyone has experience with something similar I’d love to hear ideas!
I’ve been looking into using numba as well, so if anyone has experience with that too that would be great.
I’ve not dealt much with the sparse systems in scipy but I think @cookpat4 has some experience in this realm. Coincidentally he also has experience with numba
Depending on what you want to do with numba, it could be as easy as adding the jit decorator before the function you want to speed-up, but there may be some caveats depending on how complex your code is.
I’m pretty sure at this point the issue is related to the weird-shaped convergence region or RK. I think it makes sense to just do the tried and true R-matrix variational approach; e.g. solve the Bloch-Schroedinger equation on a Lagrange mesh, rather than the Schroedinger equation itself in the coordinate basis, That way we are explicitly solving a small, dense linear system of the form Ax=b, which is much more numerically nice than a general operator equation of the form F[phi] = 0.
