Tridiagonal Systems

Charles Kankelborg

2001-01-21

1 The Problem

A matrix equation Mu = r is just a way of writing N linear, algebraic equations. The ith equation is:

∑_{j} M_{ij} u_{j} = r_{i}.  (1)

The standard problem is to solve for u given M and b. A square matrix is said to be tridiagonal when

M_{ij} = 0, |i − j| > 1.  (2)

For example,

(\begin{matrix}b_{1} & c_{1} & 0 & 0 \\ a_{2} & b_{2} & c_{2} & 0 \\ 0 & a_{3} & b_{3} & c_{3} \\ 0 & 0 & a_{4} & b_{4} \end{matrix}) (\begin{matrix}u_{1} \\ u_{2} \\ u_{3} \\ u_{4} \end{matrix}) = (\begin{matrix}r_{1} \\ r_{2} \\ r_{3} \\ r_{4} \end{matrix}) .  (3)

2 Algorithm

Have a look Numerical Recipes Section 2.6. The Numerical Recipes function tridag is easily understood without LU decomposition.

To solve the 4 × 4 system in equation 3 for u, we begin by using row 1 to eliminate a_{2} from row 2. Similarly, we eliminate all the a_{i} row by row.

Here is how the tridag implementation works. The first row stands for the equation:

b_{1} u_{1} + c_{1} u_{2} = r_{1}.

Let β_{1} = b_{1},  (4)
v_{1} = \frac{r_{1}}{β_{1}} ,  (5)
γ_{2} = \frac{c_{1}}{β_{1}}  (6)

If we divide row 1 by β_{1}, its equation becomes

u_{1} + γ_{2} u_{2} = v_{1}.  (7)

The equation for the second row is:

a_{2} u_{1} + b_{2} u_{2} + c_{2} u_{3} = r_{2}.

Let β_{2} = b_{2} − a_{2} γ_{2},
v_{2} = \frac{(r_{2} − a_{2} v_{1})}{β_{2}},
γ_{3} = \frac{c_{2}}{β_{2}} .

We now subtract a_{2} times row 1 from row 2. Dividing the result by β_{2}, our new row 2 is written:

u_{2} + γ_{3} u_{3} = v_{2}.  (8)

This process may be repeated over the range i = 2, 3, … N, using the recursion formulas

γ_{i} = \frac{c_{i−1}}{β_{i−1}} ,  (9)
β_{i} = b_{i} − a_{i} γ_{i},  (10)
v_{i} = \frac{r_{i} − a_{i} v_{i−1}}{β_{i}} .  (11)

The result is that the matrix equation is rewritten in uper triangular form. For our 4 × 4 example,

(\begin{matrix}1 & γ_{2} & 0 & 0 \\ 0 & 1 & γ_{3} & 0 \\ 0 & 0 & 1 & γ_{4} \\ 0 & 0 & 0 & 1 \end{matrix}) (\begin{matrix}u_{1} \\ u_{2} \\ u_{3} \\ u_{4} \end{matrix}) = (\begin{matrix}v_{1} \\ v_{2} \\ v_{3} \\ v_{4} \end{matrix}) .  (12)

This problem is trivially solved by starting from the bottom and back substituting as you go up:

u_{N} = v_{N},  (13)
u_{i} = v_{i} − γ_{i+1} u_{i+1}.  (14)

3 How to Break the Algorithm

If β_{i} = 0, for some i, the program quits to avoid a division by zero. The authors of Numerical Recipes assure us that this doesn’t happen very often.

It can be shown that β_{i} = 0 cannot occur in the calculation of y′′_{i} for the cubic spline.