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)

In Mathematica, to invert a matrix, one uses the command

Inverse[matrixa]

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

Click here to see example

 

 

 

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.

Click here to see example

 

 

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

 

* Read coefficients and constants

 

      PRINT *, 'ENTER NUMBER OF EQUATIONS'

      READ *, N

      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.

Click here to see example