Introduction

There are two kinds of finite element equation solvers : the direct and
iteration methods.

The direct method is based on the Gaussian elimination method , and its original
form was

modified by various ways to reduce the number of operations in computing , to
improve the

matrix condition number for accuracy, and to reduce the storage space of the
coefficient

matrix. Roughly speaking the total number of operations to finish up the
Gaussian

elimination procedure to solve a system of linear equations with n number of
equations, is

the order of n3 , that is, O(n3), while the required space to store the
coefficient matrix is

O(n2) if all of the coefficients are stored . If a coefficient is stored in 32
bit form, then we

need 64n2 bits space. In other words, if we have 20MB space for the coefficient
matrix,

then the largest n is about 790. This means that we can solve only 790
equations, and then

if a cubic hexagonal structure is modeled by 8 node slid brick elements, we can
only

decompose it into 5x5x5 meshes. This is not acceptable setting in practice. If
we have

20MB space for the coefficients, we would like to solve at least 10,000
equations. We can

find considerably many researches in the 1970s and 80s to reduce the storage
space

required. Typical outcome from such heavy research in the area of the finite
element

method and computational science, are, for example, the band method, skyline (
or profile )

method, wave front ( frontal ) method, and others. Especially, the wave front
method that

was developed by B. Iron was the highest achievement of the finite element
method in the

early 1970s. Because of this program to solve a system of linear equations that
is specially

designed by using the nature of linear elasticity, the finite element method
could attract wide

range people for its application, and we started solving “large scale” problems
with

O(1000) equations using 128K main memory computers with magnetic tape drives
which

provide leap frog type advancement in engineering analysis. He basically
developed an

algorithm which only requires storage of non- zero coefficients . It is noted that
the

coefficient matrix, that is the global stiffness matrix in the finite element
analysis, contains

mostly zeros. In other words, there are considerably less number of non-zero
terms than

zeros. Thus, as Iron developed, if we have an algorithm of the Gaussian
elimination

method which requires only non-zero terms, we can reduce large amount of storage
space

as well as computing time by skipping for zeros which do not make any influence.
After

Iron’s work, there were many modifications and variations, and then the skyline
( profile )

method became popular among researchers and students because of its simplicity.
A short

coming of the wave front method is its complexity of algorithm and programming,

although it was considerably simplified later, and it was applied by very
specialists of the

finite element method. On the other hand, the skyline method published in the
book by

Bathe was so simple that any engineering students could understand if they have

background of the Gaussian elimination method, and they could make their own
computer

programs without much effort. The skyline method requires more space than the
wave

front method, but it only requires storage of non-zero and zero terms bounded by
the socalled

skyline that is the most outer limit of non-zeros in columns of the coefficient
matrix.

It can be regarded as a variation of the band method that stores the non-zero
terms in a

rectangular matrix.

Direct Methods

The original Gaussian elimination method without pivoting
can be explained as follows.

Let a discrete finite element equation ( which is a typical matrix equation, and
represents a

system of linear equation with n number of unknowns and n number of equations )
:

Ax = b

or

The first step is elimination of the terms
of the first column using the first

row of the coefficient matrix A under the assumption that
is not equal to zero by using

the algorithm

This operation leads the new matrix equation

Here note that although we are using same characters
with the original matrix equation,

their values are different since we replaced the original ones by the algorithm
described.

The next is to eliminate the terms of the
second column by the similar

algorithm

and obtain the new matrix equation

Repeat this process until the upper triangular coefficient
matrix is obtained :

so that it can be solved from the last equation step by
step :

and

The process making the upper triangular coefficient matrix
is called the forward

elimination, the process modifying the right hand side is called the reduction,
and the last

step to solve the system linear equations with the upper triangular coefficient
matrix is

called the back substitution . If we examine this process, and if the coefficient
matrix has the

special structure such that all of non-zero terms are located within the band
from the

diagonals :