Page 61 - Numerical Methods for Chemical Engineering

50 1 Linear algebra

and
v 0 = v x (y = 0) v N+1 = v x (y = B) (1.252)

To enforce the no-slip boundary conditions, we merely set in the right-hand side vector of
(1.250)

v 0 = 0 v N+1 = V up (1.253)
By this technique, we have converted the original differential equation into a set of algebraic
equations for the values of the velocity at each grid point.
The matrix in (1.250) has a special structure: the only nonzero elements are located along
the main diagonal and on the diagonals immediately above and below. Such a matrix is said
to be tridiagonal. This special structure allows us to solve this system in a fraction of the
time required by brute-force Gaussian elimination.
Remember that the number of FLOPs required to perform Gaussian elimination for a
3
system of N equations is 2N /3. If we want, in general, to obtain an accurate solution of a
differential equation, we may need to use a grid of 100 or more points, so that the number
of FLOPs required is on the order of one million. In addition to CPU time, the memory
required to store the matrix is significant. A matrix for a system with N unknowns contains
2
N elements, each requiring its own location in memory. As N increases, these numbers
become much larger. The number of FLOPs required to perform full Gaussian elimination
on a system of 1000 unknowns is on the order of one billion, and storing the matrix requires
one million locations in memory.
Because here nearly all of the matrix elements are zero, much of the effort of brute-
force Gaussian elimination is a waste of time. As the matrix is tridiagonal, we only need
to perform one row operation per column to zero the element immediately below the diag-
onal. Also, since each row only contains three non zero values, the number of FLOPs
required for each row operation is a small number, not on the order of N as is the case
generally. Thus, we can remove two of the nested for loops of the Gaussian elimina-
tion algorithm so that the total number of FLOPs scales only linearly with N. This is a
very important point, as this trick of taking advantage of the tridiagonal structure makes
feasible the numerical solution of this system even when the number of grid points is
large.

Banded and sparse matrices

More generally, a matrix is said to be sparse if most of its elements are zero. A sparse
matrix can be stored in memory efficiently by recording only the positions and values of
its nonzero elements. Let us say that we have an N × N matrix A with at most only three
2
non-zero elements per row. While the total number of elements is N , at most only 3N
2
are non-zero, and the matrix is sparse for large N.If N nz N is the number of nonzero
elements in A, or at least an upper bound to the number of nonzero elements, we can store
all the information that we need about A in two integer vectors of length N nz (for the row
and column positions of the nonzero elements) and a single real vector of length N nz for
their values. For example, the matrix