# solve Ax=b for outrigger A matrix python

I implement Crank-Nicolson 2D finite-difference method.

I get a matrix A which is banded with 1 band above and below the main diagonal, but also contains 2 additional bands , further apart from the main diagonal, so it is NOT penta-diagonal.

A picture showing the structure is below. My matrix is the RHS one. The LHS is easy, it's the penta-diagonal one. I couldn't find up until now a way to solve Ax = b with A being the RHS matrix from the photo in python. I could barely find a name for it, in these lecture notes https://ocw.mit.edu/ans7870/2/2.086/F12/MIT2_086F12_notes_unit5.pdf it is called an 'outrigger' matrix (page 403).

At the moment I am using spsolve from from scipy.sparse.linalg, into which I feed two arguments, namely sparse.csc_matrix(A) and sparse.csc_array(b), where A and b have been defined initially as A = sparse.dok_matrix((size, size), dtype=np.complex64) and b = sparse.dok_array((size, 1), dtype=np.complex64), then populated with values by iterating element by element through them. It is extremely slow and I was wondering maybe someone more experienced knows a way to exploit the structure appearing in A.

Thank you!

• Using a sparse matrix data type for A is good, but your right-hand side vector b is probably dense, so I would just use a normal numpy array for that. There are also better ways to construct A than going through the dok format first, for example, the diags method. Feb 21 at 17:51
• Thank you! Much appreciated. Having a look at it now. Feb 21 at 18:09
• What's the size of your matrix? Feb 21 at 21:54
• These are small matrices. Every finite element solver solves these sorts of matrices in less than a minute on typical laptops today. Feb 22 at 18:09
• There's an example in github.com/mandli/numerical-methods-pdes/blob/master/… on how to use the spdiags command. This is for 2D elliptic problems but should not be hard to figure out your problem. Feb 23 at 15:46