Section 4.3 Pivoting
The elementary implementation of Gaussian Elimination in Section 4.2 suffers two main problems:- It fails if the element
A(j,j)is zero when we evaluateL(i,j). - There might be a loss of significant digits in the difference
A(i,k)-L(i,j)*A(j,k)ifL(i,j)is very big.
\begin{equation*}
A=\begin{pmatrix}
\epsilon&1\\
1&1\\
\end{pmatrix}.
\end{equation*}
After the first step, the element A(2,2) becomes \(1-1/\epsilon\text{.}\) If \(\epsilon\) is small enough, then \(1/\epsilon\) is so much larger than 1 that \(1-1/\epsilon=-1/\epsilon\,\text{,}\) namely we get the same decomposition as if the original matrix were
\begin{equation*}
\begin{pmatrix}
\epsilon&1\\
1&0\\
\end{pmatrix}.
\end{equation*}
A concrete example. Let us make concrete calculations in an example in the floating point system \(D_3\) with base 10 and three significant digits that we introduced in Section 1.2. Consider the system
\begin{equation*}
Ax=b
\end{equation*}
\begin{equation*}
\begin{pmatrix}
0.00025&1\\
1&1\\
\end{pmatrix}
\begin{pmatrix}
x\\
y\\
\end{pmatrix}
=
\begin{pmatrix}
1\\
2\\
\end{pmatrix}
\end{equation*}
Its exact solution is \(x=\frac{4000}{3999},y=\frac{3998}{3999}\text{.}\) Notice that these numbers are both equal to \(1.00\) in \(D_3\). Using the LU decomposition, though, we get
\begin{equation*}
L
=
\begin{pmatrix}
\phantom{400}1&0\\
4000&1\\
\end{pmatrix}
\end{equation*}
and
\begin{equation*}
U
=
\begin{pmatrix}
0.00025&1\\
0&1-4000\\
\end{pmatrix}
=
\begin{pmatrix}
0.00025&1\\
0&-4000\\
\end{pmatrix}
\end{equation*}
and the right-hand-side becomes
\begin{equation*}
L^{-1}b = \begin{pmatrix}\phantom{-400}1&0\\ -4000&1\\ \end{pmatrix}\begin{pmatrix}1\\ 2\\ \end{pmatrix}
=\begin{pmatrix}1\\ 2-4000\\ \end{pmatrix}=\begin{pmatrix}1\\ -4000\\ \end{pmatrix}
\end{equation*}
(recall that we are in a 3-digits system!) In other words, the system is equivalent to
\begin{equation*}
\left\{
\begin{array}{c}
0.00025x+y=1\\
-4000y=-4000\\
\end{array}
\right.
\end{equation*}
so that \(y=1.00\) and therefore \(x=0.00\text{.}\) Notice that we got a relative error of 100% in \(x\text{!!!}\)
Pivoting. This kind of problems can be fixed by rearranging at every step the rows of \(A\) so that, when we evaluate L(i,j), we always divide by the largest possible number. We illustrate this method by applying it to the system of the previous example:
\begin{equation*}
\begin{pmatrix}
0.00025&1\\
1&1\\
\end{pmatrix}.
\begin{pmatrix}
x\\
y\\
\end{pmatrix}
=
\begin{pmatrix}
1\\
2\\
\end{pmatrix}
\end{equation*}
As discussed above, dividing by \(0.00025\) is not a good idea, so we re-write it as
\begin{equation*}
\begin{pmatrix}
1&1\\
0.00025&1\\
\end{pmatrix}.
\begin{pmatrix}
x\\
y\\
\end{pmatrix}
=
\begin{pmatrix}
2\\
1\\
\end{pmatrix}
\end{equation*}
Clearly nothing has changed: the solution of a system does not depend on the order we use to write its equations! Now we use LU decomposition:
\begin{equation*}
L=
\begin{pmatrix}
1&0\\
0.00025&1\\
\end{pmatrix},\;
U=
\begin{pmatrix}
1&1\\
0&1\\
\end{pmatrix}
\end{equation*}
In other words, the system becomes
\begin{equation*}
\begin{pmatrix}
1&1\\
0&1\\
\end{pmatrix}.
\begin{pmatrix}
x\\
y\\
\end{pmatrix}
=
L^{-1}
\begin{pmatrix}
2\\
1\\
\end{pmatrix}
=
\begin{pmatrix}
2\\
1\\
\end{pmatrix}
\end{equation*}
(remember that we are in a 3 digits floating-point system!) whose solution is clearly \(x=y=1.00\text{.}\) The error we experienced before pivoting disappeared completely!
Pivoting code. Below is an implementation of the LU method with pivoting. As in the last implementation of the LU method in , we use a single for loop. Notice that, since while pivoting we need to exchange rows in P, A and L, we need to define from the beginning all matrices. We define L to be the zero matrix and P to be the unit matrix.