Use of scipy sparse in ode solver

Question

I am trying to solve a differential equation system

x´=Ax 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

How large and sparse is A?

It looks like A is just a constant in this function:

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

Is vector 1d? Then np.transpose(vector) does nothing.

For calculation purposes you want

Asparse = sparse.csr_matrix(A)

Does np.dot(Asparse, vector) work? np.dot is supposed to be sparse aware. If not, try Asparse*vector. This probably produces a dense matrix, so you may need (Asparse*vector).A1 to produce a 1d array.

But check the timings. Asparse needs to quite large and very sparse to perform faster than A in a dot product.