ME 203

• Introduction
• Connect
• Apply
• Reflect
• Extend

This module refers to files and folders from the tar file module5.tar.

Solving Linear Equations - The Gaussian Elimination Method
There are a number of existing web sites that provide excellent tutorials on the Gaussian method; try this one from PurpleMath and this definition and basic description from Wikipedia.

Let's consider the following 3 linear equations:

 4x - 2y +  z =15
-3x -  y + 4z = 8
  x -  y + 3z =13

We've solved equations like these enough times over the years to know that our tasks are:

• In order to solve for 3 unknowns, we need 3 equations.
• In order to solve for one of the variables, we need to eliminate, or zero out, the other 2 variables.
  • We do this two equations at a time, by multiplying the lower of the two equations by a number that will give the coefficient of the 1st variable to eliminate, and then subtracting the bottom equation from the top one.
  • We do this step again for the 1st and last equations, once again choosing a multiple that will give the coefficient we desire, and then subtracting.
  • We then repeat this process for the next variable, y, this time centering on (Or pivoting) around the middle equation and eliminating the variables from the top and bottom equations.
  • We do it one last time for the last variable, z, making the bottom equation the pivot.
    
  So, for the first two of our 3 equations, we multiply the 1st equation by 3, and the 2nd by 4, so that their x coefficients both equal 12. We'll then add the equations to eliminate the x variable.
  
  3( 4x - 2y +  z = 15)
  4(-3x -  y + 4z =  8) which becomes
  
    12x - 6y + 3z = 45
   -12x - 4y +16z = 32
     0x -10y +19z = 77
  

• Swapping out the 2nd equation for this one gives us this:
  
     4x - 2y +  z = 15
     0x -10y +19z = 77
      x -  y + 3z = 13
  
  to remove the x variable in the 3rd equation, we can multiply the 3rd equation by -4 and then add it to the 1st equation:
  
        4x - 2y +  z =15
    -4(  x -  y + 3z =13)  which becomes
  
     4x - 2y +  z = 15
    -4x + 4y -12z =-52
     0x + 2y -11z =-37
  
  plugging this back in as our 3rd equation gives us this:
  
     4x - 2y +  z = 15
     0x -10y +19z = 77
     0x + 2y -11z =-37
  
  
• Now that the x variable has been eliminated in the bottom 2 equations, we can work with them to eliminate the y variable. If we multiply the 3rd equation by 5, it will have a 10y, while the 2nd has a -10y, so all we'll have to do is add them together to get 0y:
  
     -10y + 19z = 77
   5(  2y - 11z =-37)    while becomes
  
     -10y + 19z =  77
      10y - 55z =-185
       0y - 36z =-108    or
  
           -36z =-108,   and dividing each side be -36 gives us
  
              z = 3
  
  
• at this point, we back substitue 3 for z in the 2nd equation to find for y:
  
    -10y + 19z   = 77
  
    -10y + 19(3) = 77
  
    -10y +  57   = 77
  
    -10y         = 77 - 57
  
    -10y         = 20
  
       y         = -2
  
  
• And finally, we plug in z = 3 and y = -2 into the 1st equation:
  
      4x - 2y + z    = 15
  
      4x - 2(-2) + 3 = 15
  
      4x + 4 + 3     = 15
  
      4x + 7         = 15
  
      4x             = 15 - 7
  
      4x             = 8
  
       x             = 2
  

Gaussian elimination differs from the above method in 2 ways:

• Instead of finding a multiple for, say, x's coefficient in equation a to make it equal to x's coefficient in equation b, we use the ratio of the 2 coefficients. This mechanism for variable elimination makes far more sense from a programmatic point of view, which is our intention, anyway.
• In order to guard against dividing by zero, the Gaussian method rearranges the equations so that larger coefficients are on the diagonal (i.e., variable x's coefficient for equation 1, variable y's coeefficient for equation 2, and variable z's coefficient in equation 3, which makes the diagonal). This mechanism of rearranging larger coefficients onto the diagonal is called pivoting.

We'll implement the Gaussian method with one change; instead of pivoting, which adds a significant amount of complexity to the program, we'll use an if-then conditional to test for, and avoid, division by zero.

Also keep in mind that in our program, gauss_no_pivot.c, we refer to a variable called pivot; this is the coefficient for the variable on the diagonal, and is not the same as pivoting.

Let's take a closer look at our set of equations:

 4x - 2y +  z = 15
-3x -  y + 4z =  8
  x -  y + 3z = 13

First off, we notice that since the variables are in the same place in each equation, we can effectively remove them without affecting the mathematics:

 4 - 2 + 1 = 15
-3 - 1 + 4 =  8
 1 - 1 + 3 = 13

In addition, as long as we conserve positive and negative signs, we can remove arithmetic symbols:

 4  -2   1    15
-3  -1   4     8
 1  -1   3    13

Finally, notice that we have what looks like an 3x3 matrix next to a 3 element vector:

As a result, we can operate on the coefficients of each equation's variables as elements within an array, using the programmatic tools the C programming language provides for such arrays.

When you decompress the module5.tar file, and go into the module5/ directory, you'll notice a subdirectory called gauss. Within this subdirectory is the program gauss_no_pivot.c that implements the Gaussian elimination method.

Let's step through the code:

• Lines 1 thru 17 are generic include files and variable declarations. The variables coefficient, pivot and ratio (Line 9) are the most interesting ones; ratio will be calculated based on coefficient and pivot to eliminate the variable for the current row, or equation (Line 70) as we discussed above.
  
  
• Lines 19 to 54 are pretty standard issue for arrays as well,
  • 19 to 22 is the the if-then to confirm existence of the files,
  • Lines 24 thru 35 make sure that matrix/vector sizes match,
  • 38 & 39 define the vector b and the matrix A based on the scanned sizes.
  • Lines 41 thru 51 scan the matrix and vector data from files and assigns it into the appropriate array locations (And prints the matrix to screen as well).
  • We close the files after the scans in lines 53 & 54.
    
    
• The doubly nested for loops in lines 59 thru 81 do most of the heavy lifting in this program:
  • The variable pivot is assigned to the coefficient on the diagonal (Line 60).
  • Instead of pivoting equations based on coefficient size to avoid division by zero, we test to see if the divisor is zero in the if-then statement of lines 62 thru 64. If it is, we exit.
  • The 2nd for loop initiated in line 67 starts on the row (jrow) directly below the current one (irow).
    • In line 70, we define ratio as the coefficient of the variable we're currently trying to zero out (Line 69) divided by the pivot from line 60.
    • The 3rd, innermost for loop initiated in line 73 is where we reset each row element (Basically, each coefficient of the equation) by subtracting from it, the ratio multiplied by the element (Line 75).
  • We also must calculate this formula for the right-hand side equality value, which we do in line 79.
    
    
• At the completion of the doubly nested for loops, the last row will contain the coefficient for one variable, and the vector at that row (n-1) will contain what it equals; if we divide both the coefficient and the right-hand side by coefficient, we find out what that one variable is equal to. Line 83 performs this division and resets the vector at the last row to this value.
  
  
• Now, we must go back up through each row, plugging in the variable values we've calculated for the rows (equations) below it, to find for all variables. This is what we do in the nested for loop in lines 88 to 96.
  • Line 88 defines the outer loop as going from the second to last row (We already found for the last row in line 83) back up to the first row.
  • The inner for loop in lines 90 through 92 takes each variable already found, multiplies it by the value of the coeffiecient for that variable, and continually removes that from the value on the right-hand side of the equality in vector b.
    
    
• At the end of the loop in line 95, we divide what we calculated in that inner loop by the coefficient at the diagonal to finish the calculation of the right-hand side value.
  
  
• Once the above nested for loop is finished, we have all values for all variables in the vector b.
  • The for loop at the end of the program in lines 98 to 100 outputs the values for each equation from vector b.

The Lower/Upper (LU) Method
For an in-depth description of the Lower/Upper method (As well as another take on the Gaussian method) check out the efunda (Engineering Fundamentals) page covering the subject. We provide a basic definition below.

The Lower/Upper method for solving n linear equations breaks the matrix of coefficients up into 2 triangles; a Lower triangle that the bottom left of the matrix, and an Upper triangle that covers the top right of the matrix.

< graphic of the Upper & Lower triangles here >

So the equation now becomes

Ax = b => LUx = b

and if we assign y = Ux, then we only have to solve the following 2 systems:

Ly = b
Ux = y

This turns out to be more efficient than the Gaussian method above; The amount of time needed to solve using Gaussian decomposition will be based on the cube of the number of equations, whereas back-substitution of a triangular system of equations will take around n-squared time.

Let's take a look at the implementation of Lower/Upper decomposition, lu.c, inside the LU directory within the module5 tar file.

• Lines 6 to 16 declare the variables of the program, including
  • incrementing counter integers in line 7,
  • 10 x 10 matrices containing the coefficients of the linear equations and their results in line 9 and 13 (The actual matrix size is dictated within the file containing the matrix data, and can be no larger than 10 due to static nature of the matrix declaration),
  • vectors of size 10 to hold the equation data and temporary scratch pad calculations during the variable computations, in lines 10 thru 12, and 14 thru 16.
• Lines 18 through 36 open the matrix.dat and vector.dat files, and read their data into matrix A and vector b, respectively.