Tridiagonal Matrices: Thomas Algorithm

1 Tridiagonal Matrices: Thomas Algorithm W. T. Lee . MS6021, Scientific Computation, University of Limerick The Thomas Algorithm is an efficient way of solving Tridiagonal matrix systems. It is based on LU decompo- sition in which the matrix system M x = r is rewritten as LU x = r where L is a lower triangular matrix and U is an upper triangular matrix . The system can be efficiently solved by setting U x = and then solving first L = r for and then U x = for x. The Thomas Algorithm consists of two steps. In Step 1 decomposing the matrix into M = LU and solving L = r are accomplished in a single downwards sweep, taking us straight from M x = r to U x = . In step 2 the equation U x = is solved for x in an upwards sweep. I. STAGE 1 Row 2 a2 x1 + b2 x2 +c2 x3 = r2.

2 A2 Row 1 a2 x1 + a2 1 x 2 = a2 1. In the first stage the matrix equation M x = r is converted New Row 2 (b2 a2 1 ) x2 +c2 x3 = r2 a2 1. to the form U x = . Initially the matrix equation looks like: Divide through by (b2 a2 1 ) to get . b1 c1 0 0 0 0 x1.. r1 c2 r2 a2 1. x1 + x2 =. a2 b2 c2 0 0 0 x2 r2 b2 a 2 1 b2 a 2 1. 0 a3 b3 c3 0 0 x3 r3 .. = We can rewrite this as 0 0 a4 b4 c4 0 x4 r4 .. 0 0 0 a5 b5 c5 x5 r5 c2 r2 a2 1. x2 + 2 x3 = 2 , 2 = , 2 = . 0 0 0 0 a6 b6 x6 r6 b2 a 2 1 b2 a 2 1. Row 1.. 1 1 0 0 0 0 x1 1. b1 x1 + c1 x2 = r1 0. 1 2 0 0 0 x2 2 .. 0 a3 b3 c3 0 0 x3 r3 .. Divide through by b1 = . 0 0 a4 b4 c4 0 x4 r4 .. c1 r1 . x1 + x2 = 0 0 0 a5 b5 c5 x5 r5 .. b1 b1. 0 0 0 0 a6 b6 x6 r6. Rewrite: Row 3. c1 r1. x1 + 1 x2 = 1 , 1 = , 1 =.

3 B1 b1 a3 x2 + b3 x3 + c3 x4 = r3. Use a3 times row 2 of the matrix to eliminate the first term 1 1 0 0 0 0 x1 1.. a2 b2 c2 0 0 0 x2 r2 a3 (x2 + 2 x3 = 2 ). 0 a3 b3 c3 0 0 x3 r3 .. = . 0 0 a4 b4 c4 0 x4 r4 .. 0 Row 3 a3 x2 + b3 x3 +c3 x4 = r3. 0 0 a5 b5 c5 x5 r5 . a3 Row 2 a3 x2 + a3 2 x 3 = a3 2. 0 0 0 0 a6 b6 x6 r6. New Row 3 (b3 a3 2 ) x3 +c3 x4 = r3 a3 2. Row 2. Divide through by (b3 a3 2 ) to get a2 x1 + b2 x2 + c2 x3 = r2 c3 r3 a3 2. x2 + x3 =. Use a2 times row 1 of the matrix to eliminate the first term b3 a 3 2 b3 a 3 2. a2 (x1 + 1 x2 = 1 ) We can rewrite this as c3 r3 a3 2. x3 + 3 x4 = 3 , 3 = , 3 = . b3 a 3 2 b3 a 3 2. ; 2.. 1 1 0 0 0 0 x1 1 1 1 0 0 0 0 x1 1. 0 1 2 0 0 0 x 0 1 2 0 0 0 x . 2 2 2 2 . 0 0 1 3 0 0 x3 3 0 0 1 3 0 0 x3 3.

4 = = . 0 0 a4 b4 c4 0 0 0 0 1 0 . x4 r4 4 x4 4 .. 0 0 0 a5 b5 c5 x5 r5 . 0 0 0 0 1 5 x5 5 .. 0 0 0 0 a6 b6 x6 r6 0 0 0 0 a6 b6 x6 r6. Row 4. Row 6. a4 x3 + b4 x4 + c4 x5 = r4 a6 x5 + b6 x6 = r6. Use a4 times row 3 of the matrix to eliminate the first term Use a6 times row 5 to eliminate the first term. a4 (x3 + 3 x4 = 3 ) a6 x5 + a6 5 x6 = a6 5. Row 4 a4 x3 + b4 x4 +c4 x5 = r4 Resulting in a4 Row 3 a4 x3 + a4 3 x 4 = a4 3. (b6 a6 5 ) x6 = r6 a6 5. New Row 4 (b4 a4 3 ) x4 +c4 x5 = r4 a4 3. Divide through by b6 a6 5 to get Divide through by (b4 a4 3 ) to get r6 a6 5. c4 r4 a4 3 x6 = 6 , 6 =. x3 + x4 = b6 a 6 5. b4 a 4 3 b4 a 4 3. We can rewrite this as . 1 1 0 0 0 0 x1 1. c4 r4 a4 3 0 1 2 0 0 0 x . x4 + 4 x5 = 4 , 4 = , 4 =.

5 2 2 . b4 a 4 3 b4 a 4 3 0.. 0 1 3 0 0 x3 3 .. = . 0 0 0 1 4 0 x4 4 .. 0 0 0 0 1 5 x5 5 . 1 1 0 0 0 0 x1 1. 0 1 2 0 0 0 x 0 0 0 0 0 1 x6 6. 2 2 . 0 0 1 3 0 0 x3 3 .. = At this point the matrix has been reduced to upper diagonal 0 0 0 1 4 0 x4 4 form, so our equations are in the form U x = .. 0 0 0 a5 b5 c5 x5 r5 .. 0 0 0 0 a6 b6 x6 r6. II. STAGE 2. Row 5. The matrix equation is now in a form which is trivial to a5 x4 + b5 x5 + c5 x6 = r5 solve for x. We start with the last row and work our way up. The final equation is already solved Use a5 times row 4 of the matrix to eliminate the first term x6 = 6 . a5 (x4 + 4 x5 = 4 ). Row 5 a5 x4 + b5 x5 +c5 x6 = r5.. 1 1 0 0 0 0 x1 1. a5 Row 4 a5 x4 + a5 4 x 5 = a5 4 0. 1 2 0 0 0 x . 2 2.

6 New Row 5 (b5 a5 4 ) x5 +c5 x6 = r5 a5 4 0.. 0 1 3 0 0 x3 3 .. = . 0 0 0 1 4 0 . Divide through by (b5 a5 4 ) to get x4 4 .. 0 0 0 0 1 5 x5 5 . c5 r5 a5 4 0 0 0 0 0 1 x6 6. x4 + x5 =. b5 a 5 4 b5 a 5 4. We can rewrite this as c5 r5 a5 4. x5 + 5 x6 = 5 , 5 = , 5 = . b5 a 5 4 b5 a 5 4. 3. Row 5: x5 + 5 x6 = 5 . Row 1: x1 + 1 x2 = 1 . Rearrange to get: x5 = 5 5 x6 . Rearrange to get: x1 = 1 1 x2 .. 1 1 0 0 0 0 x1 1 1 1 0 0 0 0 x1 1. 0 1 2 0 0 0 x 0 1 2 0 0 0 x . 2 . 2 2 . 2. 0 0 1 3 0 0 x3 3 0 0 1 3 0 0 x3 3 .. = = . 0 0 0 1 4 0 0 0 0 1 0 . x4 4 4 x4 4 .. 0 0 0 0 1 5 x5 5 0 0 0 0 1 5 x5 5 . 0 0 0 0 0 1 x6 6 0 0 0 0 0 1 x6 6. Row 4: x4 + 4 x5 = 4 . At this point x, the solution to the matrix equation, is fully Rearrange to get: x4 = 4 4 x5.

7 Determined. III. IN PRACTICE. 1 1 0 0 0 0 x1 1. 0 1 2 0 0 0 x The Thomas Algorithm is used because it is fast and be- 2 . 2. 0.. 0 1 3 0 0 x3 3 . cause Tridiagonal matrices often occur in practice. (This ar- = gument is slightly circular because people often manipulate 0 0 0 1 4 0 x4 4 .. 0 the problems they are working on to reduce them to solving a 0 0 0 1 5 x5 5 . Tridiagonal matrix problem.) Although it is rare, the Algorithm 0 0 0 0 0 1 x6 6 can be unstable if bi ai i 1 is zero or numerically zero for any i. This will occur if the Tridiagonal matrix is singular, but Row 3: x3 + 3 x4 = 3 . in rare cases can occur if it is non-singular. The condition for Rearrange to get: x3 = 3 3 x4 . the Algorithm to be stable is . 1 1 0 0 0 0 x1 1 kbi k > kai k + kci k 0 1 2 0 0 0 x2 2.

8 0.. 0 1 3 0 0 x3 3 . for all i. The matrix problems which result from the discreti- = sation of partial differential equations nearly all satisfy this 0 0 0 1 4 0 x4 4 .. 0 criterion. 0 0 0 1 5 x5 5 .. If the Algorithm is numerically unstable then you must re- 0 0 0 0 0 1 x6 6 arrange the equations: known as pivoting. Standard LU de- composition algorithms for full or banded matrices include Row 2: x2 + 2 x3 = 2 . pivoting. (But first you should check to make sure you have Rearrange to get: x2 = 2 2 x3 . not made a mistake in formulating the problem.).. 1 1 0 0 0 0 x1 1. 0 1 2 0 0 0 x2 2 .. 0 0 1 3 0 0 x3 3 .. = . 0 0 0 1 4 0 x4 4 .. 0 0 0 0 1 5 x5 5 .. 0 0 0 0 0 1 x6 6.