Use of scipy sparse in ode solver

Question

I am trying to solve a differential equation system

$$x´=Ax\quad \text{with } x(0) = f(x)$$

in Python, where $A$ indeed is a complex sparse matrix.

For now i have been solving the system using the scipy.integrate.complex_ode class in the following way.

def to_solver_function(time,vector):
    sendoff = np.dot(A, np.transpose(vector))
    return sendoff

solver = complex_ode(to_solver_function)
solver.set_initial_value(f(x),0)

solution = [f(x)]
for time in time_grid:
    next = solver.integrate(time)
    solution.append(next)

This has been working OK, but I need to "tell the solver" that my matrix is sparse. I figured out that I should use

Asparse = sparse.lil_matrix(A)

but how do I change my solver to work with this?

Answer 1

For Scipy’s ODE module, the function that you feed it (in your case to_solver_function) as a blackbox that it provides with a state and that returns a vector. It does not care about what happens inside it, in particular it never touches A.

So, the only thing that you need to do is to ensure that your matrix–vector multiplication inside to_solver_function is aware of the sparseness of your matrix. I am not very familiar with scipy.sparse, but there is a paragraph of documentation on this and it tells me that the following should do the job:

sendoff = A.dot(vector)

---

^{Sidenote and blatant self-advertising: If you want to speed up things, I wrote a wrapper around Scipy’s ODE that just-in-time-compiles the derivative (in your case to_solver_function) to speed up things. Given that in your case, the most computation time should be spent on the matrix–vector multiplication and you do that already pretty efficiently, it may give you as much of a boost. Also, it does not handle complex numbers yet, so you would have to do that manually.}