Matrices (and Simultaneous Linear Equations)

One of the most powerful aspects of computers and their applications in science is their ability to reduce the usually tedious manipulations surrounding the use of matrices to relatively simple operations.  Because of this, many scientists can, for the first time, appreciate the power and beauty of matrix operations.

The following graphic shows some simple properties of matrix operations.

These operations can easily be carried out using a computer.  For example, in FORTRAN 90, the multiplication of two matrices is a single command

C = MATMUL(A,B)

# Orvis shows the following example of matrix inversion (Fig 10-1)

The natural application of matrix techniques is in the solution of sets of linear equations.

The previous Excel example showed how to do this operation in Excel.  The following example shows how to do this operation in Mathematica using several different commands.

Now we can ask how we might do this problem in FORTRAN.  To do this, we will work an example by hand to show how a particular method works.

This procedure for solving systems of linear equations is called the Gauss-Jordan method.  It can be illustrated schematically in the following flow diagram which will serve as a template for creating the FORTRAN program to implement it.

What the Gauss-Jordan routine does is to invert the A matrix (the coefficient matrix), multiply it by the inverse matrix to get the I matrix.  This method gives both b and the inverse of A.  It is a relatively slow method.

The following graphics shows two different FORTRAN programs that implement the Gauss-Jordan method.

Program #1

Program #2

A faster procedure that does not produce an inverse matrix is the method of Gauss elimination treated in Nyhoff example 11.1

PROGRAM LINSYS

************************************************************************

* Program to solve a linear system using Gaussian elimination.         *

* Identifiers used are:                                                *

*     LIMROW : maximum number of rows in the matrix                    *

*     LIMCOL : maximum number of colums (LIMROW + 1) in the matrix     *

*     N      : number of equations and unknowns                        *

*     I, J   : subscripts                                              *

*     LIN    : matrix for the linear system                            *

*     X      : solution                                                *

*     SINGUL : indicates if system is (nearly) singular                *

*                                                                      *

* Input:   The number of equations, the coefficients, and the          *

*          constants of the linear system                              *

* Output:  The solution of the linear system or a message indicating   *

*          that the system is (nearly) singular                        *

************************************************************************

INTEGER LIMROW, LIMCOL

PARAMETER (LIMROW = 10, LIMCOL = LIMROW + 1)

DOUBLE PRECISION LIN(LIMROW, LIMCOL), X(LIMROW)

INTEGER N, I, J

LOGICAL SINGUL

PRINT *, 'ENTER NUMBER OF EQUATIONS'

DO 10 I = 1, N

PRINT *, 'ENTER COEFFICIENTS AND CONSTANT OF EQUATION ',

+            I, ': '

READ *, (LIN(I,J), J = 1, N + 1)

10    CONTINUE

* Use subroutine GAUSS to find the solution,

* and then display the solution

CALL GAUSS(LIN, LIMROW, LIMCOL, N, X, SINGUL)

IF (.NOT. SINGUL) THEN

PRINT *, 'SOLUTION IS'

DO 20 I = 1, N

PRINT 100, I, X(I)

100         FORMAT(1X, 'X(', I2, ') =', F8.3)

20       CONTINUE

ELSE

PRINT *, 'MATRIX IS (NEARLY) SINGULAR'

END IF

END

**GAUSS*****************************************************************

* Subroutine to find solution of a linear system of N equations in N   *

* unknowns using Gaussian elimination, provided a unique solution      *

* exists.  The coefficients and constants of the linear system are     *

* stored in the matrix LIN, which has LIMROW rows and LIMCOL columns.  *

* If the system is singular, SINGUL is returned as true, and the       *

* solution X is undefined.  Local identifiers used are:                *                        *     I,J,K  : subscripts                                              *

*     MULT   : multiplier used to eliminate an unknown                 *

*     ASBPIV : absolute value of pivot element                         *

*     PIVROW : row containing pivot element                            *

*     EPSIL  : a small positive real value ("almost zero")             *

*     TEMP   : used to interchange rows of matrix                      *

*                                                                      *

* Accepts: Two-dimensional array LIM, integers LIMROW, LIMCOL, and N   *

* Returns: One-dimensional array X and logical value SINGUL            *

************************************************************************

SUBROUTINE GAUSS(LIN, LIMROW, LIMCOL, N, X, SINGUL)

DOUBLE PRECISION LIN(LIMROW, LIMCOL), X(LIMROW), TEMP, MULT, EPSIL

PARAMETER (EPSIL = 1D-15)

INTEGER N, PIVROW

LOGICAL SINGUL

SINGUL = .FALSE.

DO 50 I = 1, N

*        Locate pivot element

ABSPIV = ABS(LIN(I,I))

PIVROW = I

DO 10 K = I + 1, N

IF (ABS(LIN(K,I)) .GT. ABSPIV) THEN

ABSPIV = ABS(LIN(K,I))

PIVROW = K

END IF

10       CONTINUE

*        Check if matrix is (nearly) singular

IF (ABSPIV .LT. EPSIL) THEN

SINGUL = .TRUE.

RETURN

END IF

*        It isn't, so interchange rows PIVROW and I if necessary

IF (PIVROW .NE. I) THEN

DO 20 J = 1, N + 1

TEMP = LIN(I,J)

LIN(I,J) = LIN(PIVROW,J)

LIN(PIVROW,J) = TEMP

20          CONTINUE

END IF

*        Eliminate Ith unknown from equations I + 1, ..., N

DO 40 J = I + 1, N

MULT = -LIN(J,I) / LIN(I,I)

DO 30 K = I, N + 1

LIN(J,K) = LIN(J,K) +  MULT * LIN(I,K)

30          CONTINUE

40       CONTINUE

50    CONTINUE

* Find the solutions by back substitution

X(N) = LIN(N, N + 1) / LIN(N,N)

DO 70 J = N - 1, 1, -1

X(J) = LIN(J, N + 1)

DO 60 K = J + 1, N

X(J) = X(J) - LIN(J,K) * X(K)

60       CONTINUE

X(J) = X(J) / LIN(J,J)

70    CONTINUE

END

Another technique that can be used for large matrices is the Gauss-Seidel method, which is described in the following graphic.

The Gauss-Seidel method is illustrated in Orvis, chapter 10 by the following example.