Let $\mathbf{U}^n = [u^n_0, u^n_2, ..., u^n_{M-1}]^T.$ The Crank Nicolson method gives the following equation with two matrices: $$\mathrm{A}\mathbf{U}^{n+1} = \mathrm{B}\mathbf{U}^n.$$

$\mathrm{A}$ consists of diagonals $-\alpha$, $(1+2\alpha)$, and $-\alpha$, from upper to lower, respectively. $\mathrm{B}$ consists of $\alpha$, $(1-2\alpha)$, and $\alpha$, again from upper to lower diagonals respectively. Here $\alpha = D\frac{\Delta t}{2(\Delta x)^2}$. More explicitly:

$$ \mathrm{A} = \begin{bmatrix} a & b & \\ -\alpha & (1+2\alpha) & -\alpha \\ & -\alpha & (1+2\alpha) & -\alpha \\ & & \ddots & \ddots & \ddots\\ & & & -\alpha & (1+2\alpha) & -\alpha\\ & & & & c & d\\ \end{bmatrix} $$

$$ \mathrm{B} = \begin{bmatrix} f & g & \\ \alpha & (1-2\alpha) & \alpha \\ & \alpha & (1-2\alpha) & \alpha \\ & & \ddots & \ddots & \ddots\\ & & & \alpha & (1-2\alpha) & \alpha\\ & & & & h & j\\ \end{bmatrix} $$

Where the top and bottom rows, $a,b,c,d,f,g,h,j$ will have to be adjusted for boundary conditions.

Ref [https://math.stackexchange.com/questions/1259983/how-do-i-turn-this-crank-nicolson-type-equation-into-three-vectors-representing]

Question: It seems like this is normally solved by Gaussian elimination. Given that the matrix $A$ remains the same, however, would it not be more prudent to simply determine $$\mathbf{U}^{n+1} = \mathrm{A}^{-1}\mathrm{B}\mathbf{U}^n \equiv \mathrm{C}\mathbf{U}^n$$ and then multiply by $C$ for every step? The initial upfront cost of determining $C=A^{-1}B$ might be large, but if this is to be run for thousands of steps it seems like it might be more efficient. Or perhaps the resulting matrix $A^{-1}C$ is ugly enough (i.e. few zeros) that it is more efficient to just do Gaussian elimination on the left-hand-side?


1 Answer 1


Your priority should be to always use software that is sufficiently accurate, easy to maintain and sufficiently fast for the task at hand.

You are computing approximations to the solution of a differential equation. Your largest source of error is likely to be the discretization error obtained by replacing the differential equation with a system of linear equations. In order to retain the ability to estimate the discretization error using Richardson extrapolation, you need to solve the linear systems with an error that is somewhat smaller than the discretization error. However, since the discretization error is normally much larger than the rounding error, you need not worry too much about how you solve the linear system. In particular, when given the choice between recycling an LU factorization or using the $C = A^{-1} B$ it is unlikely that accuracy will ever be a concern.

Computing $C = A^{-1}B$ can be accomplished by computing an LU factorization of $A$ and solving the linear system $AC=B$ with respect to $C$. Hence there is no real difference between maintaining a code that does repeated multiplication with $C$ and a code that recycles the LU factorization again and again.

In LAPACK you have routines that can exploit the banded structure of your linear system. The function dgttrf will compute the LU factorization of a general tridiagonal matrix using Gaussian elimination with partial pivoting. The runtime is $O(n)$. The function dgttrs will use the factorization computed by dgttrf to solve a linear system. The runtime is $O(nm)$ where $m$ is the number of columns of the right-hand side matrix.

We can now consider the issue of speed. Explicitly forming the matrix $C = A^{-1} B$ has one significant downside compared with recycling the LU factorization of $A$. You extend your storage requirement from $O(n)$ to $O(n^2)$. Does it matter? Yes. You need to store the vector $U^n$, the matrix B (three diagonals), the LU factorization of A (three diagonals) and the output vector $U^{n+1}$. This this is roughly 8n words of memory. If n=2^9 = 512 and you are using double precision floating point arithmetic, then the cost is $64n = 32$ kB memory, a common size for the L1 data cache of many computers. If you explicitly form $C$ as an $n$-by-$n$ matrix, then you need $8n^2 = 2^21$ = 2 MB of memory just to store the matrix. You will have to hold the matrix in the L3 cache of your computer. The latency of the L1 cache is significantly smaller than the latency of the L3 cache, hence explicitly forming $C =A^{-1}B$ is likely to slow you down significantly.

So, to write code that it accurate, easy to maintain and not unnecessarily slow, you will use library routines to compute a banded LU factorization of $A$ and you will recycle the LU factorization again and again.

If this is still too slow, then you will have to abandon any hope of having software that is easy to maintain, but this is a question for another day.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .