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. enter image description here

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!

  • 4
    $\begingroup$ 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. $\endgroup$ Feb 21, 2022 at 17:51
  • $\begingroup$ Thank you! Much appreciated. Having a look at it now. $\endgroup$
    – velenos14
    Feb 21, 2022 at 18:09
  • $\begingroup$ What's the size of your matrix? $\endgroup$ Feb 21, 2022 at 21:54
  • 2
    $\begingroup$ These are small matrices. Every finite element solver solves these sorts of matrices in less than a minute on typical laptops today. $\endgroup$ Feb 22, 2022 at 18:09
  • 1
    $\begingroup$ 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. $\endgroup$ Feb 23, 2022 at 15:46


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.