5
$\begingroup$

Below there are two algorithms for solving tridiagonal linear systems of the form $$ \left[ \begin{array}{ccccc|c} b_1 & c_1 & & & &d_1\\ a_2 & b_2 & c_2 & & & d_2\\ & \ddots & \ddots & \ddots & & \vdots\\ & & a_{n-1} & b_{n-1} & c_{n-1} & d_{n-1}\\ & & & a_n & b_n & d_n \end{array} \right]. $$ I called them Algorithms A and B. Both of them are equivalent to Gaussian elimination, but with important difference in the form of the resulting triangular (bidiagonal) matrix.

My main question is: which one of them is more preferrable?

Algorithm A is the one that described in Wikipedia and many textbooks, it is called Thomas algorithm and is implemented, for example, in Numerical Recipes in some tricky form. Algorithm B is more straightforward and, in my opinion, is more numerically stable in cases when $|b_i|\gg|a_i|+|c_i|$ . Though I haven't seen Algorithm B in texbooks, note that exactly this algorithm is implemented in the mentioned Wikipedia article, see "Implementation in Fortran 90", while "Implementation in Matlab" deals with Algorithm A ("Implementation in C" in its current state is a mess that does not seem to work at all).


$$ \begin{array}{|c|c|}\hline \mathbf{Algorithm\ A} & \mathbf{Algorithm\ B}\\\hline \textit{% Elimination}&\textit{% Elimination}\\ \begin{array}{l} \tilde c_1=c_1/b_1\\ \tilde d_1=d_1/b_1\\ \mathbf{for }\quad i=2 \quad \mathbf{to}\quad n-1 \quad \textbf{do}\\ \quad q=b_i-a_i c_{i-1}\\ \quad \tilde c_i=c_i/q\\ \quad \tilde d_i=(d_i-a_i \tilde d_{i-1})/q\\ \mathbf{end do}\\ \tilde d_n=(d_n-a_n \tilde d_{n-1})/(b_{n}-a_n \tilde c_{n-1})\\ \\ \end{array} & \begin{array}{l} \\ \\ \hat b_1=b_1\\ \mathbf{for }\quad i=2 \quad \mathbf{to}\quad n \quad \textbf{do}\\ \quad q=a_i/\hat b_{i-1}\\ \quad \hat b_i=b_i-q c_{i-1}\\ \quad \hat d_i=d_i-q \hat d_{i-1}\\ \mathbf{end do}\\ \\ \\ \end{array}\\ \hline \textit{% Resulting system} & \textit{% Resulting system}\\ \left[ \begin{array}{ccccc|c} 1 & \tilde c_1 & & & &\tilde d_1\\ & 1 & \tilde c_2 & & & \tilde d_2\\ & & \ddots & \ddots & & \vdots\\ & & & 1 & \tilde c_{n-1} & \tilde d_{n-1}\\ & & & & 1 & \tilde d_n \end{array} \right] & \left[ \begin{array}{ccccc|c} \hat b_1 & c_1 & & & &\hat d_1\\ & \hat b_2 & c_2 & & & \hat d_2\\ & & \ddots & \ddots & & \vdots\\ & & & \hat b_{n-1} & c_{n-1} & \hat d_{n-1}\\ & & & & \hat b_n & \hat d_n \end{array} \right]\\ \hline \textit{% Backsubtitution} & \textit{% Backsubtitution}\\ \begin{array}{l} \\ x_n=\tilde d_n\\ \mathbf{for }\quad i=n-1 \quad \mathbf{downto}\quad 1 \quad \textbf{do}\\ \quad x_i=\tilde d_i-\tilde c_i x_{i+1}\\ \\ \end{array} & \begin{array}{l} \\ x_n=\hat d_n/\hat b_n\\ \mathbf{for }\quad i=n-1 \quad \mathbf{downto}\quad 1 \quad \textbf{do}\\ \quad x_i=(\hat d_i-c_i x_{i+1})/\hat b_i\\ \\ \end{array}\\\hline \end{array} $$

$\endgroup$
  • $\begingroup$ I looked around a bit at some of the primary resources online, and I couldn't find any papers in the obvious places (LAPACK Working Notes, etc...) discussing this particular routine, and the LAPACK routine for this, xgtsl, still bears Jack's original copyright. Neither of your approaches employs pivoting, which is probably a more important factor than the other differences between them. $\endgroup$ – Aron Ahmadia Apr 19 '12 at 12:05
  • $\begingroup$ @AronAhmadia: That's true, unless the tridiagonal system is diagonally dominant, in which case no partial pivoting is necessary. $\endgroup$ – Paul Apr 20 '12 at 4:44
3
$\begingroup$

Both algorithms compute $LU$ decompositions (solving against $L$ while it is being formed) and then solve against the resulting $U$. The difference is that Algorithm A forces $U$ to have a diagonal of all ones (we say that $U$ is unit-diagonal), while Algorithm B forces $L$ to have a unit diagonal (this is the usual convention).

Regardless of whether or not one is more stable than the other, both are a bad idea; as @AronAhmadia mentioned, you should use an algorithm which performs partial pivoting. I would go with the LAPACK routine dgtsv.

$\endgroup$
  • $\begingroup$ Thank you! Your answer is quite surprising for me, since I've always believed Numerical Recipes (apps.nrbook.com/rollover/index.html, section 2.4) which says that "tridiagonal algorithm is the rare case of the algorithm that, in practice, is more robust than theory says it should be". Please look at two more questions in the updated post. $\endgroup$ – faleichik Apr 20 '12 at 11:25
  • $\begingroup$ @faleichik - I'm going to roll back your edit, as it's rather unfair to answerers to modify a question to ask new things (the point of modifying a question is to improve clarity, not get more things). $\endgroup$ – Aron Ahmadia Apr 20 '12 at 13:58
  • $\begingroup$ If you want to ask a new question, that's fine, though you should just look up the definition of pivoting, there are plenty of easily-constructed cases where you create a divide-by-zero or divide-by-near-zero problem. Regarding the article describing the implementation, there may not be much, since it is a fairly trivial extension of LU. $\endgroup$ – Aron Ahmadia Apr 20 '12 at 14:00
  • $\begingroup$ Honestly, I don't think that my updating was dishonest. Firstly I wanted to add the questions you've deleted as a comments to Jack's answer, but decided to add them to the main post in order to attract more attention from other people. I don't really think that these additional questions deserve a separate post. $\endgroup$ – faleichik Apr 20 '12 at 14:11
  • $\begingroup$ As for the pivoting: you know, I'm aware of how it works for tridiagonal matrices, the comments in 'dgtsv' are pretty clear. I was just wondering if someone considered this problem in a textbook because all texts I've seen deal with Thomas algorithm and nothing more. $\endgroup$ – faleichik Apr 20 '12 at 14:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.