25.3.5.7 Substitution

We introduced the LU-decomposition as a method for solving a linear system of equations. However, until now we have only decomposed the matrix A into L and U. Obviously, we now want to take the second step and see that this matrix decomposition helps us in solving a linear system of equations.

For this we return to the system of Eq. 25.2 from which we also took the matrix. The right-hand side is

$\overrightarrow{b}=\left(\begin{array}{c}2\\ 4\\ 0\end{array}\right)$

We will now use the function DoLUDecomposition in combination with Eq. 25.20 and Eq. 25.21 for solving the linear system of equations

$\begin{array}{c}\mathbf{L}\overrightarrow{y}=\overrightarrow{b}\\ \mathbf{U}\overrightarrow{x}=\overrightarrow{y}\end{array}$

described by Eq. 25.18 and Eq. 25.19.

Calculating $\overrightarrow{y}$. We first solve Eq. 25.19 by applying Eq. 25.20 for finding $\overrightarrow{y}$ by forward-substitution. For this, we find the following equations

$\begin{array}{ll}{y}_1=b_1& =2\\ y_2=b_2-l_{21}{y}_1& =4-\frac{7}{3}\cdot 2=-\frac{2}{3}\\ y_3=b_3-{\displaystyle \sum _{k=1}^2{l}_{3kyk}=-\left(l_{31y1}{+}_{32y2}\right)=-\frac{1}{3}\cdot 2+\frac{1}{19}\cdot \frac{2}{3}}& =-\frac{12}{19}\end{array}$

Calculating $\overrightarrow{x}$. Now that we have $\overrightarrow{y}$ we proceed by finding $\overrightarrow{x}$ by applying Eq. 25.21. For this we find the following equations

$\begin{array}{lllll}{x}_3& =& \frac{1}{u_{33}}\left(y_3\right)=-\frac{19}{10}\frac{12}{19}& =& -\frac{6}{5}\\ x_2& =& \frac{1}{u_{22}}\left(y_2-u_{23}x3\right)=-\frac{3}{19}\left(-\frac{2}{3}-\frac{11}{3}\cdot \frac{6}{5}\right)& =& \frac{4}{5}\\ x_1& =& \frac{1}{u_{11}}\left(y_1-u_{12}{x}_2-u_{13}{x}_3\right)=\frac{1}{3}\left(2-4\cdot \frac{4}{5}-2\cdot \frac{6}{5}\right)& =& -\frac{6}{5}\end{array}$

which obviously is the correct result (see Eq. 25.5). As you can see, we are able to find the solution to our original system of equation pretty much effortlessly once we have performed the LU-decomposition. In this case, we did not even have to perform a row permutation and there seems to be no problem with stability. This is obviously not the case, as our solver would still fail if we found a zero on the downwards-diagonal due to the division in the back-substitution performed in Eq. 25.21. As we have already done previously, this problem can be circumvented by changing the order of the equations, i.e., by implementing pivoting.

Maple Worksheet. Obviously, at this point we want to implement this scheme as a function. The function DoLUSubstitution implements the substitution steps required for solving linear systems of equations using LU-decomposition. It requires the matrices L and U as well as the vector $\overrightarrow{b}$ as parameters. It consists of two loops, the first of which performs the forward-substitution (see line 6) while the second one performs the back-substitution (see line 14).

For the forward-substitution, the sum required for Eq. 25.20 is created in line 9 before applying Eq. 25.20 in line 11 to find the individual entries of $\overrightarrow{y}$. For the back-substitution, we first calculate the sum required for Eq. 25.21 in line 17 before applying Eq. 25.21 in line 19 for finding $\overrightarrow{x}$.

If we run our example (see line 24) we obtain the correct result

$\overrightarrow{x}=\left(\begin{array}{c}0.4\\ 0.8\\ -1.2\end{array}\right)$

Stability Consideration. As we can see, the LU-decomposition gives us the correct result and seems to be sufficiently stable. However, this is only because the example we chose happens to have no zeros on the downwards-diagonal. As already discussed before, these would result in a division by zero. As sparsely populated matrices with a lot of zero entries are very common, it is a good idea to implement pivoting for LU-decomposition as well.

25.3.5.8 Pivoting

As we have already done for SOR in section 25.3.4.4, we will implement partial pivoting for LU-decomposition. This method is often referred to as permutating LU-decomposition (PLU). For this we change the implementation of the function DoLUDecomposition as shown in listing 25.7. As you can see, we use a pivoting vector piv that is first initialized to the sequence (1, 2, . . . , n) in line 8. Inside of the line loop line 12 we perform the pivoting operation as before (see line 14). We keep the original order for the matrices L and U, but we then use the pivoting vector to access the respective row in m as can be seen in line 26 and line 34. We then need to return the pivoting vector in addition to the matrices L and U in line 37 because we require the pivoting vector to correctly associate the individual entries of the vectors$\overrightarrow{x}$ and $\overrightarrow{b}$ when solving the linear system of equations.

 1 DoLUSubstitution:=proc(L::Matrix,U::Matrix,b::Vector) local n,i,k,x,y,summed:
 2   n:=Size(b,1):
 3   y:=Vector(1..n):
 4   x:=Vector(1..n):
 5   #forward-substitution - find y[i]
 6   for i from 1 to n do
 7     summed:=0:
 8     for k from 1 to 1-1 do
 9       summed:=summed+L[i,k]*y[k]:
10     end do:
11     y[i]:=evalf(b[i]-summed):
12   end do:
13   #back-substitution - find x[i]
14   for i from n by -1 to 1 do
15     summed:=0:
16     for k from i+1 to n do
17       summed:=summed+U[i,k]*x[k]:
18     end do:
19     x[i]:=evalf(1/U[i,i]*(y[i]-summed)):
20   end do:
21   return x:
22 end proc:
23
24 mml:=Matrix([[3,4,2],[7,3,1],[1,1,1]]);b:=Vector([2,4,0]);
25 DoLUSubstitution(DoLUDecomposition(m),b);

 1 DoLUDecompositionPivot:=proc(mml::Matrix) local n,i,j,k,l,summed,L,U,piv,store:
 2   n:=Size(b,1):
 3   L:=Matrix(n,n):
 4   for i to n do
 5     L[i,i]:=1:
 6   end do:
 7   U:=Matrix(n,n):
 8   piv:=Vector(1..n):
 9   for i to n do
10     piv[i]:=i:
11   end do:
12   for i to n do
13    for l from i+1 to n do
14     if abs(m[piv[l],i])>abs(m[piv[i],i]) then
15       store:=piv[l]:
16       piv[l]:=piv[i]:
17       piv[i]:=store:
18     end if:
19    end do:
20    #calculate l[ij]
21    for j from 1 to i-1 do
22     summed:=0:
23     for k from 1 to j-1 do
24       summed:=summed+U[k,j]*L[i,k]:
25     end do:
26     L[i,j]:=1/U[j,j]*(m[piv[i],j]-summed):
27    end do:
28    #calculate u[ij]
29    for j from i to n do
30     summed:=0:
31     for k from 1 to i-1 do
32       summed:=summed+U[k,j]*L[i,k]:
33     end do:
34     U[i,j]:=m[piv[i],j]-summed:
35    end do:
36   end do:
37   return piv,L,U;
38 end proc:
39
40 mml:=Matrix([[3,4,2],[7,3,1],[1,1,1]]);
41 DoLUDecompositionPivot(m);
42
43 mml:=Matrix([[0,4,2],[7,0,1],[1,1,0]]);
44 DoLUDecompositionPivot(m);

Example. If we apply the function DoLUDecompositionPivot to our original matrix (see line 40) you can see from the output that the algorithm has performed a pivoting operation by switching the first and the second rows. We have found previously that pivoting is not necessary for this problem. However, switching the two rows places the 7 on the downwards-diagonal in the first row. This makes our algorithm more stable even though it is not required.

A Second Example. In the next step, we modify our original problem by setting (pretty much deliberately) two zeros on the downwards-diagonal. If we ran this second example (see line 43) with the original implementation of our LU-decomposition using the function DoLUDecomposition, we would receive an error due to the division by zero. If we run this example with the LU-decomposition using pivoting, we avoid this error and obtain the solution

$\mathbf{L}=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ \frac{1}{7}& \frac{1}{4}& 1\end{array}\right),\mathbf{U}=\left(\begin{array}{ccc}7& 0& 1\\ 0& 4& 2\\ 0& 0& -\frac{2}{14}\end{array}\right)$

The pivoting vector we also receive shows the sequence (2, 1, 3) indicating that the algorithm switched the lines such that the zeros are moved off the downwards-diagonal.

Adapting the Function DoLUSubstitution. Now that we have adapted the LU-decomposition to support pivoting, we also need to adapt the function DoLUSubstitution. The implementation of the method DoLUSubstitutionPivot is shown in listing 25.8 and is in analogy to the implementation of LUSubstitution (see listing 25.6). The only difference is that the function expects the pivoting vector piv as additional parameter. We then use this vector for accessing the correct entries in $\overrightarrow{b}$ (see line 11). We can now run our original example (see line 24) which would obviously return the correct result again. In addition we can now also solve the second example we applied when testing the LU-decomposition with pivoting in listing 25.7. This is the output of line 27. The solution to this matrix is

 1 DoLUSubstitutionPivot:=proc(piv::Vector,L::Matrix,U::Matrix,b::Vector) local n,i,k,x,y,summed:
 2   n:=Size(m,1):
 3   y:=Vector(1..n):
 4   x:=Vector(1..n):
 5   #forward-substitution - find y[i]
 6   for i from 1 to n do
 7     summed:=0:
 8     for k from 1 to i-1 do
 9       summed:=summed+L[i,k]*y[k]:
10     end do:
11     y[i]:=evalf(b[piv[i]]-summed):
12   end do:
13   #back-substitution - find x[i]
14   for i from n by -1 to 1 do
15     summed:=0:
16     for k from i+1 to n do
17       summed:=summed+U[i,k]*x[k]:
18     end do:
19     x[i]:=evalf(1/U[i,i]*(y[i]-summed)):
20   end do:
21   return x:
22 end proc:
23
24 mml:=Matrix([[3,4,2],[7,3,1],[1,1,1]]);b:=Vector([2,4,0]);
25 DoLUSubstitutionPivot(DoLUDecompositionPivot(m),b);
26
27 mml:=Matrix([[0,4,2],[7,0,1],[1,1,0]]);
28 DoLUDecompositionPivot(m);

$\overrightarrow{x}=\left(\begin{array}{c}0.333\\ -0.333\\ 1.666\end{array}\right)=\left(\begin{array}{c}\frac{1}{2}\\ -\frac{1}{3}\\ \frac{5}{3}\end{array}\right)$

which is the correct result.

25.3.5.9 Summary

In this section we have discussed LU-decomposition, a very useful method for solving systems of linear equations. Besides being an exact method that does not require iteration, LU-decomposition allows us to use the decomposed matrices L and U for various values of $\overrightarrow{b}$ , and thus to “reuse” a solution to linear systems with changing right-hand sides. We have seen that LU-decomposition by itself is already reasonably stable. By implementing partial pivoting, we can use this method to solve almost any kind of linear system of equations, which makes it a very potent and widely applicable method.

25.3.6 Sequential Dense LU-Decomposition

25.3.6.1 Introduction

In this section we will briefly introduce a modification of the classical LU-decomposition which is a little more memory efficient. In many cases, it is not necessary to have both matrices U and L separately. As they are both triangular matrices, they can be stored in a single matrix without losing information. In fact, there is a modification to the classical LU-decomposition that works with the matrix A directly transferring it to a matrix ${\mathbf{A}}^{\text{'}}$ of the form

${\mathbf{A}}^{\text{'}}=\mathbf{L}-{\mathit{I}}_{\text{n}}+\mathbf{U}$

with I_n being the identity matrix of rank n. This scheme is referred to as sequential dense LU-decomposition. It is a very memory-efficient scheme and therefore lends itself well to implementation for numerical solvers.

25.3.6.2 Derivation

We will derive the scheme using the generic matrix A that we gradually transform to the dense LU-decomposed matrix ${\mathbf{A}}^{\text{'}}$.

$\mathbf{A}=\left(\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{3n}& a_{3n}& \cdots & a_{nn}\end{array}\right)=\left(\begin{array}{cccc}1& 0& \cdots & 0\\ l_{21}& 1& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ l_{n1}& l_{n2}& \cdots & 1\end{array}\right)\cdot \left(\begin{array}{cccc}{u}_{11}& u_{12}& \cdots & u_{1n}\\ 0& u_{22}& \cdots & u_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & u_{nn}\end{array}\right)=\mathbf{L}\cdot \mathbf{U}$

First Row. We now change the entries of A row-by-row figuring out how we need to modify them in order to end up with ${\mathbf{A}}^{\text{'}}$. We begin with the first row for which we need to find the values of u₁₁ to u₁₂. For this we find

$\begin{array}{ccc}{u}_{11}& =& a_{11}\\ u_{12}& =& a_{12}\\ u_{ij}& =& a_{1j}\\ \cdots & & \\ u_{1n}& =& a_{1n}\end{array}$

which means that the values of a₁₁ to a_1n can remain unchanged as they already represent the correct values of u₁₁ to u_1n. In algorithmic notation we could note down that

$u_{1j}=a_{1j}$

for i = 1 as the row index and i ≤ j ≤ n as the column index.

Remainder of First Column. After having cleared the first row, we proceed the rest of the first column for which we find the following equations

$\begin{array}{ccccccc}{l}_{21}{u}_{11}& =& a_{21}\to l_{21}{a}_{11}& =& a_{21}\to l_{21}a& =& \frac{a_{21}}{a_{11}}\\ l_{21}{u}_{11}& =& a_{i1}\to l_{i1}{a}_{11}& =& a_{i1}\to l_{i1}a& =& \frac{a_{i1}}{a_{11}}\\ \cdots & & & & & & \\ l_{n1}{u}_{11}& =& a_{n1}\to l_{n1}{a}_{11}& =& a_{n1}\to l_{n1}& =& \frac{a_{n1}}{a_{11}}\end{array}$

As you can see, we find the remaining values of the first column by simply dividing the entries of A by the entry in the first row and first column a₁₁. In algorithmic notation we can note

$l_{i1}=\frac{a_{i1}}{a_{11}}$

with the row index i with 1 ≤ i ≤ n. As you can see, we simply have to divide the entries a_i1 by the value a₁₁ to find l_i1. On the other hand, we may just as well simply perform this division on the matrix A itself, which will be the first modification we make. We therefore replace the first row entries of A with the values a⁰i₁ according to

$a^{\text{'}}{}_{il}=\frac{a_{il}}{a_{11}}=l_{i1}$

Second Row. At this point, the first row and the first column have been determined and the respective entries in L and U are known. We now proceed with the second row. Obviously, it would be convenient to simply repeat the process. However, we will see that this is not possible directly as the values of the entries we already determined will interfere. This is where the second modification will be made to A. By a small trick, we will change A in order to make all influences of already determined entries “disappear.” For the second row we find the following equations

$\begin{array}{ccccc}{l}_{21}{u}_{12}+u_{22}& =& a_{22}\to u_{22}& =& a_{22}-l_{21}{u}_{12}\\ l_{21}{u}_{13}+u_{23}& =& a_{23}\to u_{23}& =& a_{23}-l_{21}{u}_{13}\\ l_{21}{u}_{1j}+u_{2j}& =& a_{2j}\to u_{2j}& =& a_{2j}-l_{21}{u}_{1j}\\ \cdots & & & & \\ l_{21}{u}_{1n}+u_{2n}& =& a_{2n}\to u_{2n}& =& a_{2n}-l_{21}{u}_{1n}\end{array}$

As you can see, the operation is almost analogous to the treatment of the first row with the one exception that we cannot simply copy the values of a_2j. We need to subtract the produce l₂₁ u_1j first. Now, we could just as well simply subtract this product from a_2j and end up with the same result. In algorithmic notation, we could note down that

$u_{2j}=a_{2j}-l_{21}{u}_{1j}$

$\begin{array}{lll}{a}^{\text{'}}& =& a_{2j}-l_{21}{u}_{1j}\\ & =& a_{2j}-a^{\text{'}}{}_{21}{a}_{1j}=u_{2j}\end{array}$

where we have used Eq. 25.36 to set u_1j = a_1j and Eq. 25.37 to set l₂₁ = a⁰₂₁. Obviously, we expect a similar influence of the already determined values on all other lines as well. We will therefore perform a likewise subtraction on the other lines using the same formula according to

$\begin{array}{lll}{a}^{\text{'}}{}_{ij}& =& a_{2j}-l_{21}{u}_{1j}\\ & =& a_{ij}-a^{\text{'}}{}_{i1}{a}_{1j}=u_{ij}\end{array}$

for i ≥ 2 and j ≥ 2. We will see in a moment that we actually required this modification.

Remainder of Second Column. Having processed the second row, we now turn to the second column starting from the second row. Proceeding likewise we find the following equations

$\begin{array}{lllll}{l}_{31}{u}_{12}+l_{32}{u}_{22}& =& a_{32}\to l_{32}& =& \frac{1}{u_{22}}\left(a_{32}-l_{31}{u}_{12}\right)\\ l_{41}{u}_{12}+l_{42}{u}_{22}& =& a_{42}\to l_{42}& =& \frac{1}{u_{22}}\left(a_{42}-l_{41}{u}_{12}\right)\\ l_{i1}{u}_{12}+l_{i2}{u}_{22}& =& a_{i2}\to l_{i2}& =& \frac{1}{u_{22}}\left(a_{i2}-l_{i1}{u}_{12}\right)\\ \cdots & & & & \\ l_{n1}{u}_{12}+l_{in2}{u}_{22}& =& a_{n2}\to l_{n2}& =& \frac{1}{u_{22}}\left(a_{n2}-l_{n1}{u}_{12}\right)\end{array}$

In algorithmic notation we note

$l_{i2}=\frac{1}{u_{22}}\left(a_{i2}-i_{i1}{u}_{12}\right)$

This expression can be simplified by noting that for j = 2 Eq. 25.39 tells us that a⁰22 = u22 in which case we can rewrite the expression to

$l_{i2}=\frac{1}{a^{\text{'}}{}_{22}}\left(a_{i2}-l_{i1}{u}_{12}\right)$

As we presumed just a moment ago, we actually find that the modification we made by Eq. 25.40 was necessary and that this modification allows us to write Eq. 25.41 as

$l_{i2}=\frac{1}{a^{\text{'}}{}_{22}}\left(a_{i2}-l_{i1}{u}_{12}\right)=\frac{a_{i2}{}^{\text{'}}}{a_{22}{}^{\text{'}}}$

for j = 2. If we compare this with Eq. 25.37, we find that this is the same operation and we can change a⁰i₂ again in-place.

Iteration. Obviously, at this point we have found an iteration pattern if we repeat the process consisting of the following steps starting with A using the line index i

1. In line i replace for all (i, j) with i < j ≤: n replace $a_{ij}^i=\frac{a_{ij}}{a_{ii}}$ according to Eq. 25.37

2. Correct submatrix entries (k, j) with i < k ≤ n and i < j ≤ n: replace $a_{kj}^i=a_{kj}-a_{k,i}{a}_{i,j}$ according to Eq. 25.40

3. Repeat for next row i + 1

25.3.6.3 Example

We will quickly carry out this example using the matrix we have been using as example so far

$\mathbf{A}=\left(\begin{array}{ccc}3& 4& 2\\ 7& 3& 1\\ 1& 1& 1\end{array}\right)$

We proceed by the outlined scheme by first dividing all entries for i < j ≤ n = 3 in the first column by a₁₁ = 3 in which case we end up with

${\mathbf{A}}^{\prime }=\left(\begin{array}{ccc}3& 4& 2\\ \frac{7}{3}& 3& 1\\ \frac{1}{3}& 1& 1\end{array}\right)$

In the next step we correct the submatrix entries (k, j) for i < k ≤ n and i < j ≤ n finding

${\mathbf{A}}^{\prime }=\left(\begin{array}{ccc}3& 4& 2\\ \frac{7}{3}& -\frac{19}{3}& -\frac{11}{3}\\ \frac{1}{3}& -\frac{1}{3}& \frac{1}{3}\end{array}\right)$

We then proceed to i = 2 again dividing all values i < j ≤ n = 3 (which is only the value a₃₂) by $a_{22}=-\frac{19}{3}$ yielding

${\mathbf{A}}^{\prime }=\left(\begin{array}{ccc}3& 4& 2\\ \frac{7}{3}& -\frac{19}{3}& -\frac{11}{3}\\ \frac{1}{3}& \frac{1}{19}& \frac{1}{3}\end{array}\right)$

We then proceed by correcting the remaining submatrix which only consists of the entry a₃₃, which we correct to the final result

${\mathbf{A}}^{\prime }=\left(\begin{array}{ccc}3& 4& 2\\ \frac{7}{3}& -\frac{19}{3}& -\frac{11}{3}\\ \frac{1}{3}& \frac{1}{19}& \frac{10}{3}\end{array}\right)$

If you compare this result with Eq. 25.33, you can see that we have obtained the matrix A′ according to Eq. 25.35.

25.3.6.4 Maple Worksheet

In the next step, we will transfer the algorithm to a Maple worksheet as we will need this specific algorithm again. The listing is shown in listing 25.9. The sequential dense LU-decomposition is implemented by the function DoLUDecompositionSequentialDense. It expects the matrix to decompose as parameter m. It will modify this matrix in place, therefore it may be required to copy the matrix if the original is to be retained. First the rank of the matrix is determined (see line 2) before the individual rows are processed indexed with the variable i. For each line the columns are iterated using the variable j dividing each entry according to Eq. 25.37. After this, the submatrix is corrected. For this a new row index k is introduced. The column index j can be reused here, but the value i is required which in turn requires the new index. Line 9 implements Eq. 25.40. The last two lines (see line 16) output the example which we derived manually.

 1 DoLUDecompositionSequentialDense:=proc(mml::Matrix) local n,i,j,k:
 2   n:=Size(m,1):
 3   for i from 1 to n do
 4     for j from i+1 to n do
 5      m[j,i]:=m[j,i]/m[i,i]:
 6     end do:
 7    for k from i+1 to n do
 8     for j from i+1 to n do
 9      m[k,j]:=m[k,j]-m[k,i]*m[i,j]:
10     end do:
11    end do:
12   end do:
13   return m;
14 end proc:
15
16 mml:=Matrix([[3,4,2],[7,3,1],[1,1,1]]);
17 DoLUDecompositionSequentialDense(m);

25.3.6.5 Pivoting

As already discussed during the discussion of the regular LU-decomposition, it is generally required to use pivoting for LU-decomposition. We will therefore reimplement DoLUDecompositionSequentialDense in a way that it supports pivoting. The listing is shown in listing 25.10. We implement partial pivoting in analogy to listing 25.7 using a pivoting vector piv (see line 4). You will see that this implementation is identical to the implementation of listing 25.7. The only difference is that each access to a row i must now be made via the pivoting vector. As we only have the matrix m in which we construct the matrices L and U, we cannot keep the original order as we did for the matrices in listing 25.7.

 1 DoLUDecompositionSequentialDensePivot:=proc(mml::Matrix) local n,i,j,k,l,piv,store:
 2   n:=Size(m,1):
 3   piv:=Vector(1..n):
 4   for i from 1 to n do
 5    piv[i]:=i:
 6   end do:
 7   for i from 1 to n do
 8    for l from i to n do
 9     if abs(m[piv[l],i])>abs(m[piv[i],i]) then
10       store:=piv[l]:
11       piv[l]:=piv[i]:
12       piv[i]:=store:
13     end if:
14    end do:
15    for j from i+1 to n do
16     m[piv[j],i]:=m[piv[j],i]/m[piv[i],i]:
17    end do:
18    for k from i+1 to n do
19     for j from i+1 to n do
20      m[piv[k],j]:=m[piv[k],j]-m[piv[k],i]*m[piv[i],j]:
21     end do:
22    end do:
23   end do:
24   return piv,m;
25 end proc:
26
27 mml:=Matrix([[0,4,2],[7,3,1],[1,1,0]]);
28 DoLUDecompositionPivot(m);
29 DoLUDecompositionSequentialDensePivot(m);

If we run our example (see line 27), we will see that we obtain the correct solution

${\mathbf{A}}^{\prime }=\left(\begin{array}{ccc}0& 4& 2\\ 7& 0& 1\\ \frac{1}{7}& \frac{1}{4}& -\frac{9}{14}\end{array}\right)$

which is the same result as obtained from DoLUDecompositionPivot, keeping in mind that L and U are combined according to Eq. 25.35. Also, by inspecting piv you can see that the algorithm changed the order of the rows ending up with the sequence (2, 1, 3). The function DoLUDecompositionPivot will rearrange L and U in order to “hide” this reordering from the user. The function DoLUDecompositionSequentialDensePivot cannot do this, as the function only has a single matrix to operate on. Therefore the substitution function for DoLUDecompositionSequentialDensePivot would need to account for this inability to reorder by accessing L and U (which are then stored in m) via the pivoting vector as index.

25.3.8 Tridiagonal Matrix Algorithmml: Thomas Algorithm

25.3.8.1 Introduction

There is one algorithm which we want to introduce at this point that is commonly used to solve linear systems of equations where the matrix A is tridiagonal. This algorithm is generally referred to as the tridiagonal matrix algorithm or Thomas¹ algorithm. Please recall that tridiagonal matrices are matrices which have nonzero entries on the downwards-diagonal and in the columns right and left of the diagonal only. All other entries are zero. These matrices occur surprisingly often and they must usually not be solved by a complete LU-decomposition if the matrix is set up such that the largest values are collected on the downwards-diagonal. Again, such matrices occur frequently. We will come across one example when studying one-dimensional FEM (see Eq. 32.48). We will use this example here while introducing the algorithm. Please refer to section 32.5 for the derivation of the matrix.

25.3.8.2 Derivation

The general linear system of equation can be written as

$\begin{array}{cc}\mathbf{A}\overrightarrow{x}& =\overrightarrow{b}\\ \left(\begin{array}{cccccccc}{\beta }_1& {\gamma }_1& 0& 0& 0& \dots & 0& 0\\ {\alpha }_2& {\beta }_2& {\gamma }_2& 0& 0& \dots & 0& 0\\ 0& {\gamma }_3& {\beta }_3& {\gamma }_3& 0& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& 0\\ 0& 0& 0& 0& 0& \dots & \alpha n& {\beta }_n\end{array}\right)\cdot \left(\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\\ ⋮\\ x\left(n-1\right)\\ x_n\end{array}\right)& =\left(\begin{array}{c}{b}_1\\ b_2\\ b_3\\ b_4\\ b_5\\ ⋮\\ b\left(n-1\right)\\ b_n\end{array}\right)\end{array}$

which we may rewrite to

$\begin{array}{llll}{\beta }_1{x}_1& +{\gamma }_1{x}_2& =b_1& i=1\\ \alpha i{x}_{i-1}& +{\beta }_i{x}_i+{\gamma }_i{x}_i+1& =b_i& i=\mathrm{2..}n\end{array}$

for the i^th equation. We split up the solving process in two steps: a forward step and a back-substitution.

Forward Step: First Row. The first step is dividing the first row by β₁. Doing this, we obtain a 1 on the downwards-diagonal entry of the matrix. In this case we can rewrite Eq. 25.44 for i = 1 to

$\begin{array}{ll}{x}_1+\frac{{\gamma }_1}{{\beta }_1}{x}_2& =\frac{b_1}{{\beta }_1}\\ x_1+{{\gamma }^{\prime }}_1{x}_2& ={b^{\prime }}_2\end{array}$

The first two equations are now given by

$\begin{array}{ll}{x}_1+{{\gamma }^{\prime }}_1{x}_2& ={b^{\prime }}_1\\ {\alpha }_2{x}_1+{\beta }_2{x}_2+{\gamma }_2{x}_3& =b_2\end{array}$

Forward Step: Subsequent Rows. We now multiply the first row by α₂ and subtract it from the second row finding

$\begin{array}{ll}{x}_1+{{\gamma }^{\prime }}_1{x}_2& ={b^{\prime }}_1\\ \left({\beta }_2-{\alpha }_2{{\gamma }^{\prime }}_1\right)x_2+{\gamma }_2{x}_3& =b_2-{b^{\prime }}_1{\alpha }_2\end{array}$

We then divide the second row by (β₂ − α₂γ′₁) thus obtaining again a 1 on the diagonal entry of the second row. We then obtain

$\begin{array}{ll}{x}_1+{{\gamma }^{\prime }}_1{x}_2& ={b^{\prime }}_1\\ \hfill +x_2+\frac{{\gamma }_2}{{\beta }_2-{\alpha }_2{{\gamma }^{\prime }}_1}{x}_3& =\frac{b_2-{b^{\prime }}_1{\alpha }_2}{{\beta }_2-{\alpha }_2{{\gamma }^{\prime }}_1}\end{array}$

which we simplify to

$\begin{array}{ll}{x}_1+{{\gamma }^{\prime }}_1{x}_2& ={b^{\prime }}_1\\ \hfill +x_2+{{\gamma }^{\prime }}_2{x}_3& ={b^{\prime }}_2\end{array}$

We then proceed with all subsequent rows accordingly. The general scheme is

$\begin{array}{l}{{\alpha }^{\prime }}_i=0\\ {{\beta }^{\prime }}_i=0\\ {{\gamma }^{\prime }}_i=\frac{{\gamma }_i}{{\beta }_i-{\alpha }_i{{\gamma }^{\prime }}_{i-1}}\end{array}$

${b^{\prime }}_i=\frac{b_i-{\alpha }_i{b^{\prime }}_{i-1}}{{\beta }_i-{\alpha }_i{{\gamma }^{\prime }}_{i-1}}$

After proceeding with these modifications Eq. 25.43 can be written as

$\begin{array}{cc}\mathbf{A}\overrightarrow{x}& =\overrightarrow{b}\\ \left(\begin{array}{cccccccc}1& {{\gamma }^{\prime }}_1& 0& 0& 0& \dots & 0& 0\\ 0& 1& {{\gamma }^{\prime }}_2& 0& 0& \dots & 0& 0\\ 0& 0& 1& {{\gamma }^{\prime }}_3& 0& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& 0\\ 0& 0& 0& 0& 0& \dots & 0& 1\end{array}\right)\cdot \left(\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\\ ⋮\\ x_{(n-1)}\\ x_n\end{array}\right)& =\left(\begin{array}{c}{b^{\prime }}_1\\ {b^{\prime }}_2\\ {b^{\prime }}_3\\ {b^{\prime }}_4\\ {b^{\prime }}_5\\ ⋮\\ {b^{\prime }}_{(n-1)}\\ {b^{\prime }}_n\end{array}\right)\end{array}$

Back-Substitution. The back-substitution uses the fact that the last line of Eq. 25.48 gives the solution for x_n directly, which is simply b⁰_n. We can then substitute all the rows above the last row one by one noting that

$x_i={b^{\prime }}_i-{{\gamma }^{\prime }}_i{x}_{i+1}$

25.3.8.3 Maple Worksheet

We will quickly implement the Thomas algorithm in Maple. The function we will implement is called DoThomas and is shown in listing 25.12 together with the example from Eq. 32.48. We first restart the Maple server and include Core.txt (see line 1). The function DoThomas accepts the matrix A as parameter m and the right-hand side vector $\overrightarrow{b}$. It first determines the rank of the matrix and prepares the solution vector x (see line 4 and line 5). It then calculates the values for the first row applying Eq. 25.45 for γ′₁ and b′₁ (see line 6 and line 7). It then loops all subsequent rows applying Eq. 25.46 for finding γ′ⁱ and Eq. 25.47 for finding b′_i (see line 9 and line 10). The loop only extends to n − 1 because for the last row only b′_n can be calculated as this row has no entry for γ (see line 12). Please note that in all cases the algorithm does not set the values for α′_i and β′₁, as these are not required to find the solution.

 1 restart:read "Core.txt":
 2
 3 DoThomas:=proc(m::Matrix,b) local n,x,i:
 4   n:=Size(m,1):
 5   x:=Vector(1..n):
 6   m[1,2]:=m[1,2]/m[1,1]:
 7   b[1]:=b[1]/m[1,1]:
 8   for i from 2 to n-1 do
 9    m[i,i+1]:=m[i,i+1]/(m[i,i]-m[i,i-1]*m[i-1,i]):
10    b[i]:=(b[i]-m[i,i-1]*b[i-1])/(m[i,i]-m[i,i-1]*m[i-1,i]):
11   end do:
12   b[n]:=(b[n]-m[n,n-1]*b[n-1])/(m[n,n]-m[n,n-1]*m[n-1,i]):
13   x[n]:=evalf(b[n]):
14   for i from n-1 to 1 by -1 do
15    x[i]:=evalf(b[i]-m[i,i+1]*x[i+1]):
16   end do:
17   return x;
18 end proc:
19
20 A:=Matrix([[2 ,-1,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0],
21           [-1,2 ,-1,0 ,0 ,0 ,0 ,0 ,0 ,0],
22           [ 0,-1,2 ,-1,0 ,0 ,0 ,0 ,0 ,0],
23           [ 0, 0,-1, 2,-1,0 ,0 ,0 ,0 ,0],
24           [ 0, 0, 0,-1,2 ,-1,0 ,0 ,0 ,0],
25           [ 0, 0, 0, 0,-1,2 ,-1,0 ,0 ,0],
26           [ 0, 0, 0, 0,0 ,-1,2 ,-1,0 ,0],
27           [ 0, 0, 0, 0,0 ,0 ,-1,2 ,-1,0],
28           [ 0, 0, 0, 0,0 ,0 ,0 ,-1,2 ,-1],
29           [ 0, 0, 0, 0,0 ,0 ,0 , 0,-1,2]]):
30 b:=Vector([4,4,4,4,4,4,4,4,4,4]):
31 DoThomas(A,b);

Having completed the forward step, the algorithm then proceeds with the back-substitution first setting the value x_n directly (see line 13). It then loops (see line 14) from the second-to-last value all the way to the first setting in the values in the solution vector according to Eq. 25.49.

Having implemented DoThomas we then proceed by solving the example taken from Eq. 32.48 using the indicated values for A and $\overrightarrow{b}$ (see line 30). Executing line 31 yields the solution

${\overrightarrow{v}}_x=\left(\begin{array}{c}20\\ 36\\ 48\\ 56\\ 60\\ 60\\ 56\\ 48\\ 36\\ 20\end{array}\right)$

which is correct. In section 32.5 we will detail the use of this solution and its physical interpretation.

25.3.8.4 Numerical Stability

As with all numerical methods, it is useful to consider the criteria for numerical stability. The critical step for the numerical stability is the division operation in Eq. 25.46 and Eq. 25.47. For all cases in which the denominator results in very small values, the algorithm is prone to becoming numerically unstable. Without going into details, we will obtain a stable scheme whenever

$\left|{\beta }_i\right|\ge \left|{\alpha }_i\right|+\left|{\gamma }_i\right|$

In these cases, the denominator will be large and therefore the division operation will result in a continuous decrease in magnitude. In fact, many applications yield tridiagonal matrices which fulfill these conditions, including our example. Our implementation does not check for numerical stability. In cases where stability is a problem, LU-decomposition with pivoting should be applied which will yield stable solutions.

25.3.8.5 Summary

The Thomas algorithm is a very fast and efficient algorithm for solving linear systems of equations involving tridiagonal matrices. It is often a significantly faster alternative to using a full LU-decomposition. The Thomas algorithm is numerically quite stable and in many applications it is a good first choice to solve a tridiagonal matrix problem.