2
$\begingroup$

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!

Improve this question
$\endgroup$
6
  • 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$
    – Daniel Shapero
    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$
    – Wolfgang Bangerth
    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$
    – Wolfgang Bangerth
    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$
    – Kyle Mandli
    Feb 23, 2022 at 15:46

0

Reset to default

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.