# Ways to solve $Ax=b$ for a sparse (banded) $A$ with updates

I want to solve the time-dependent Schrodinger Equation using the Crank-Nicolson scheme. I end up with the following matrix equation

A_fixed = ...  # the matrix that does the finite difference method
B_fixed = ...  # the matrix that does the finite difference method
A = A_fixed - diags(V)  # Add the potential term V to the diagonal of A_fixed
B = B_fixed + diags(V)  # Add the potential term V to the diagonal of A_fixed

# solve for psi_new at each time step:
A.dot(psi_new) = B.dot(psi_old)   # psi_old is the solution at the previous time-step


And $A$ and $B$ are sparse matrices with all the terms along the main/off-diagonals (is this called 'banded'?). What are the best ways to solve for the above matrix equation?

I am using the scipy library sp.sparse.linalg, where it provides a solve method which uses LU decomposition. If the potential term $V$ is time-independent, then both $A$ and $B$ are actually fixed throughout the whole computation, so I can use sp.sparse.linalg.factorized to pre-factorize $A$. This gives me like a 2x speed up.

But what if $V$ is time-varying, so $A$ and $B$ changes at each time step? Are there are algorithms which work for this specific problem (where both $A$ and $B$ consist of fixed terms and some time-varying terms which are only added to the main diagonal)?

If the only non-zero entries of $A_{ij}$ have $j$ in $\{i - 1, i, i + 1\}$, then $A$ is a banded matrix with bandwidth 1. More generally, you can talk about matrices of bandwidth $k$ where $k$ is any integer. For example, if you were using a higher-order finite difference discretization that used more points to calculate a derivative, you'd get higher-bandwidth matrices.
It's possible to solve a system with a banded matrix in $O(\text{bandwidth}^2\cdot n)$ time, where $n$ is the size of the matrix; this is essentially the Thomas algorithm. Since you're using scipy already, it has a built-in banded matrix solver, although you'll probably have to change how you store the matrices to the banded format. A banded solver should be substantially faster than a general sparse factorization because it avoids having to store the factored form of the matrix.
• If I'm not mistaken, all nonzero elements in a 'banded matrix' must be 'grouped' together along the main diagonal, right? Does it mean that it only works for 1D problem, but not 2D? Since the matrix A and B are used to implement spatial derivative (finite difference), if it is 1d, then all terms will lie along/right next to the main diagonal. But if it is 2d, then when I reshape the 2d grip to a 1d vector, the finite difference along one direction must lie far away from the main diagonal, so the diagonals are 'separated'. – Physicist Aug 23 '18 at 8:29