592 Appendix A
We can now solve for thexistraightforwardly: first, we see that its last elementxN
is equal tob′′N/a′′NN. Then the remainingxiare found as
xi=
b′′i −∑N
j=i+ 1a′′ijxj
/a′′ii. (A.113)The latter procedure is calledback-substitution. It requires onlyO(N^2 )steps since
findingxirequiresO(N−i)steps.
For large matrices, the procedure described above is often unstable. This is
because the elementa′mm, which is used to zero the elements below it, may become
small. As the factor by which the lower rows are multiplied contains 1/a′mm,itis
clear that the elements of these lower rows may become very large. Combinations
of small and large numbers often cause numerical instabilities in the computer. To
avoid such problems in matrix calculations it is advised to usepivoting. Pivoting
means that the largest element (in absolute value)a′im,i=m,...,Nis looked for.
If this largest element is nota′mm, then rowsiandmare interchanged so that the
largest possible number now occurs in the denominator of the multiplication factor.
This largest element is called thepivot.
Up to now, we have described a stable and efficient method for solving the
fundamental problem of linear algebra,Eq. (A.110), findingxforonegiven vector
b. However, we may often have to face the problem of solving for several, or even
many, vectorsbi. For the case of severalbs (by several we mean less than the
dimension of the matrixA), we can perform the row changes needed to diagonalise
Aon the corresponding elements of all thebisimultaneously. It is even possible,
by initially putting thebiequal to theNunit vectors, to calculate the inverse ofA.
However, instead of following the procedure described above, it is preferable to zero
not only the part below but also above the diagonal ofA(meanwhile performing
the same changes to the elements of thebis on the right hand side) and finally to
normalise the remaining diagonal elements ofAto 1, so that it is now in unit matrix
form. Thebiare then the columns of the inverse ofA, and no back-substitution is
necessary.
For manybi(i.e. more than the dimension ofA) it is tempting to use the inverse
A−^1 ofAto solve for the solutionsxby
x=A−^1 b (A.114)but this procedure is quite susceptible to round-off errors and therefore not recom-
mended. A stable and efficient method for solvingEq. (A.110)for many vectorsbi
is theLUdecomposition. In this method, one decomposesAinto two matricesL
andUwithA=L·U. The matrixLis lower, andUupper triangular. BecauseL
andUcontain togetherN(N+ 1 )nonzero elements andAonlyN^2 , the problem is