186 6 Linear Algebra
wherea,b,care one-dimensional arrays of length1 :n. In the example of Eq. (6.15) the arrays
aandcare equal, namelyai=ci=− 1 /h^2. We can rewrite Eq. (6.15) as
Au=
b 1 c 1 0 ... ... ...
a 2 b 2 c 2 ... ... ...
a 3 b 3 c 3 ... ...
... ... ... ... ...
an− 2 bn− 1 cn− 1
an− 1 bn
u 1
u 2
...
...
...
un
=
f 1
f 2
...
...
...
fn
.
A tridiagonal matrix is a special form of banded matrix whereall the elements are zero except
for those on and immediately above and below the leading diagonal. The above tridiagonal
system can be written as
aiui− 1 +biui+ciui+ 1 =fi,
fori= 1 , 2 ,...,n. We see thatu− 1 andun+ 1 are not required and we can seta 1 =cn= 0. In many
applications the matrix is symmetric and we haveai=ci. The algorithm for solving this set of
equations is rather simple and requires two steps only, a forward substitution and a backward
substitution. These steps are also common to the algorithmsbased on Gaussian elimination
that we discussed previously. However, due to its simplicity, the number of floating point
operations is in this case proportional withO(n)while Gaussian elimination requires 2 n^3 / 3 +
O(n^2 )floating point operations. In case your system of equations leads to a tridiagonal matrix,
it is clearly an overkill to employ Gaussian elimination or the standard LU decomposition. You
will encounter several applications involving tridiagonal matrices in our discussion of partial
differential equations in chapter 10.
Our algorithm starts with forward substitution with a loop over of the elementsiand can be
expressed via the following piece of code taken from the Numerical Recipe text of Teukolsky
et al[36]
btemp = b[1];
u[1] = f[1]/btemp;
for(i=2 ; i <= n ; i++){
temp[i] = c[i-1]/btemp;
btemp = b[i]-a[i]*temp[i];
u[i] = (f[i] - a[i]*u[i-1])/btemp;
}Note that you should avoid cases withb 1 = 0. If that is the case, you should rewrite the
equations as a set of ordern− 1 withu 2 eliminated. Finally we perform the backsubstitution
leading to the following code
for(i=n-1 ; i >= 1 ; i--){
u[i] -= temp[i+1]*u[i+1];
}Note that our sums start withi= 1 and that one should avoid cases withb 1 = 0. If that is the
case, you should rewrite the equations as a set of ordern− 1 withu 2 eliminated. However, a
tridiagonal matrix problem is not a guarantee that we can finda solution. The matrixAwhich
rephrases a second derivative in a discretized form