Python: How to solve DAE with Jacobian efficiently?

Question

I am trying to use the Assimulo package to solve a set of differential algebraic equations (DAEs). I am trying to use an algorithm similar to that shown here. However, there does not seem to be an option to pass in a sparse matrix. My Jacobian matrix is very large, approximately 3000 x 3000. Do you know if there is a way to solve my DAEs more computationally efficiently?

Answer 1

In my experience with sparse ODE systems (more precisely with systems of semi-discretized PDEs), using an iterative linear solver greatly enhances numerical efficiency. As far as I know, Assimulo doesn't allow to provide a jacobian sparsity pattern, but changing the linear solver is another way to tackle this.

You would do something like:

model = Explicit_Problem(ode_function, y0=y_init, t0=t_init)
simulator = CVode(model)
sim.linear_solver = 'SPGMR'

I'm not sure if this also applies for DAE systems, but I think it's worth giving this a try.