Previous: 9.4 Arrays of Pointers
Up: 9.1 Two Dimensional Arrays
Next: 9.6 Common Errors
Previous Page: 9.4 Arrays of Pointers
Next Page: 9.6 Common Errors

9.5 An Example: Linear Algebraic Equations

As our final example program using two dimensional arrays in this chapter, we develop a program to solve systems of simultaneous linear equations. A set of linear algebraic equations, also called simultaneous equations, occur in a variety of mathematical applications in science, engineering, economics, and social sciences. Examples include: electronic circuit analysis, econometric analysis, structural analysis, etc. In the most general case, the number of equations, , may be different from the number of unknowns, ; thus, it may not be possible to find a unique solution. However, if equals , there is a good chance of finding a unique solution for the unknowns.

Our next task is to solve a set of linear algebraic equations, assuming that the number of equations equals the number of unknowns:

LINEQNS: Read the coefficients and the right hand side values for a set of linear equations; solve the equations for the unknowns.

The solution of a set of linear equations is fairly complex. We will first review the process of solution and then develop the algorithm in small parts. As we develop parts of the algorithm, we will implement these parts as functions. The driver will just read the coefficients, call on a function to solve the equations, and call a function to print the solution.

Let us start with an example of a set of three simultaneous equations in three unknowns: , , and .

We can use arrays to represent such a set of equations; a two dimensional array to store the coefficients, a one dimensional array to store the values of the unknowns when solved, and another one dimensional array to store the values on the right hand side. Later, we will include the right hand side values as an additional column in the two dimensional array of coefficients. Each row of the two dimensional array stores the coefficients of one of the equations. Since the array index in C starts at 0, we will assume the unknowns are the elements x[0], x[1], and x[2]. Similarly, the elements in row zero are the coefficients in the equation number 0, the elements in row one are for equation number one, and so forth.

Then using arrays, a general set of linear algebraic equations with unknowns may be expressed as shown below:

a[0][0] * x[0] +  a[0][1]*x[1] + ...+  a[0][m - 1] * x[m - 1]  =  y[0]
  a[1][0] * x[0] +  a[1][1]*x[1] + ...+  a[1][m - 1] * x[m - 1]  =  y[1]
  ...
  a[n-1][0]*x[0] +   ...              +  a[n-1][m - 1]*x[m - 1]  =  y[n-1]
The unknowns and the right hand side are assumed to be elements of one dimensional arrays: x[0], x[1],, x[m - 1] and y[0], y[1],, y[n - 1], respectively. The coefficients are assumed to be elements of a two dimensional array: a[i][j] for and . The coefficients of each equation correspond to a row of the array. For our discussion in this section, we assume that equals . With this assumption, it is possible to find a unique solution of these equations unless the equations are linearly dependent, i.e. some equations are linear combinations of others.

A common method for solving such equations is called the Gaussian elimination method. The method eliminates (i.e. makes zero) all coefficients below the main diagonal of the two dimensional array. It does so by adding multiples of some equations to others in a systematic way. The elimination makes the array of new coefficients have an upper triangular form since the lower triangular coefficients are all zero.

The modified equivalent set of equations in unknowns in the upper triangular form have the appearance shown below:

a[0][0]*x[0]+  a[0][1]*x[1] +    ...         +   a[0][n-1]  *x[n-1] =   y[0]
                 a[1][1]*x[1] +    ...         +   a[1][n-1]  *x[n-1] =   y[1]
                                 a[2][2]*x[2]..+   a[2][n-1]  *x[n-1] =   y[2]
                                                   a[n-1][n-1]*x[n-1]=  y[n-1]
The upper triangular equations can be solved by back substitution. Back substitution first solves the last equation which has only one unknown, x[n-1]. It is easily solved for this value --- x[n-1] = y[n-1]/a[n-1][n-1]. The next to the last equation may then be solved - since x[n-1] has been determined already, this value is substituted in the equation, and this equation has again only one unknown, x[n-2]. The unknown, x[n-2], is solved for, and the process continues backward to the next higher equation. At each stage, the values of the unknowns solved for in the previous equations are substituted in the new equation leaving only one unknown. In this manner, each equation has only one unknown which is easily solved for.

Let us take a simple example to see how the process works. For the equations:

1 * x[0] + 2 * x[1] + 3 * x[2] = 6
     2 * x[0] + 3 * x[1] + 1 * x[2] = 6
     1 * x[0] + 0 * x[1] + 2 * x[2] = 3
We first reduce to zero the coefficients in the first column below the main diagonal (i.e. array index zero). If the first equation is multiplied by -2 and added to the second equation, the coefficient in the second row and first column will be zero:

1 * x[0] + 2 * x[1] + 3 * x[2] =  6
     0 * x[0] - 1 * x[1] - 5 * x[2] = -6
     1 * x[0] + 0 * x[1] + 2 * x[2] =  3
Similarly, if the first equation is multiplied by -1 and added to the third equation, the coefficient in the third row and first column will be zero:
1 * x[0] + 2 * x[1] + 3 * x[2] =  6
     0 * x[0] - 1 * x[1] - 5 * x[2] = -6
     0 * x[0] - 2 * x[1] - 1 * x[2] = -3
Coefficients in the first column below the main diagonal are now all zero, so we do the same for the second column. In this case, the second equation is multiplied by a multiplier and added to equations below the second; thus, multiplying the second equation by -2 and adding to the third makes the coefficient in the second column zero:
1 * x[0] + 2 * x[1] + 3 * x[2] = 6
     0 * x[0] - 1 * x[1] - 5 * x[2] = -6
     0 * x[0] + 0 * x[1] + 9 * x[2] = 9
We now have equivalent equations with an upper triangular form for the non-zero coefficients. The equations can be solved backwards - the last equation gives us x[2] = 1. Substituting the value of x[2] in the next to the last equation and solving for x[1] gives us x[1] = 1. Finally, substituting x[2] and x[1] in the first equation gives us x[0] = 1.

From the above discussion, we can see that a general algorithm involves two steps: modify the coefficients of the equations to an upper triangular form, and solve the equations by back substitution.

Let us first consider the process of modifying the equations to an upper triangular form. Since only the coefficients and the right hand side values are involved in the computations that modify the equations to upper triangular form, we can work with these items stored in an array with rows and columns (the extra column contains the right hand side values).

Let us assume the process has already reduced to zero the first columns below the main diagonal, storing the modified new values of the elements in the same elements of the array. Now, it is time to reduce the lower column to zero (by lower column, we mean the part of the column below the main diagonal). The situation is shown in below:

a[0][0]   a[0][1]   ...         a[0][k]      ...          a[0][n]
     0         a[1][1]   ...         a[1][k]      ...          a[1][n]
     0         0         a[2][2]...  a[2][k]      ...          a[2][n]
     ...       ...       ...         ...          ...          ...
     0         0         0...        a[k][k]      a[k][k+1]..  a[k][n]
     0         0         0...        a[k+1][k]    ...          a[k+1][n]
     0         0         0...        a[k+2][k]    ...          a[k+2][n]
     ...       ...       ...         ...          ...          ...
     0         0         0...        a[n-1][k]    ...          a[n-1][n]
The column represents the right hand side values with a[i][n] equal to y[i]. We multiply the row by an appropriate multiplier and add it to each row with index greater than . Assuming that a[k][k] is non-zero, the row multiplier for addition to the row () is:
-a[i][k] / a[k][k]
The row multiplied by the above multiplier and added to the row will make the new a[i][k] zero. The following loop will reduce to zero the lower column:
/*   Algorithm: process_column
          Reduces lower column k to zero.
     */
     for (i = k + 1; i < n; i++) {      /* process rows k+1 to n-1 */
       m = - a[i][k] / a[k][k]          /* multiplier for ith row */

for (j = k; j <= n; j++) /* 0 thru k-1 cols. are zero */ a[i][j] += m * a[k][j]; /* add kth row times m */ /* to ith row. */ }

However, before we can use the above loop to reduce the lower column to zero, we must make sure that a[k][k] is non-zero. If the current a[k][k] is zero, all we need to do is exchange this row with any higher indexed row with a non-zero element in the column. After the exchange of the two rows, the new a[k][k] will be non-zero. The above loop is then used to reduce the lower column to zero. The non-zero element, a[k][k] used in the multiplier is called a pivot.

So, there are two steps involved in modifying the equations to upper triangular form: for each row find a pivot, and reduce the corresponding lower column to zero. If a non-zero pivot element is not found, then one or more equations are linear combinations of others, the equations are called linearly dependent, and they cannot be solved.

Figures 9.22 and 9.23 show the set of functions that convert the first rows and columns of an array to an upper triangular form. These and other functions use a user defined type, status, with possible values ERROR returned if there is an error, and OK returned otherwise. The type status is defined as follows:

typedef enum {ERROR, OK} status;
We also assume a maximum of MAX equations, so the two dimensional array must have MAX rows and MAX+1 columns. Figure 9.21 includes the header file with the defines and function prototypes used in the program.

Since precision is important in these computations, we have used formal parameters of type double. The two dimensional arrays can store coefficients for a maximum of MAX equations (rows) and have MAX + 1 columns to accommodate the right hand side values.

The function uptriangle() transforms coefficients of the equations to an upper triangular form. For each k from 0 through n-1, it calls findpivot() to find the pivot in the column. If no pivot is found, findpivot() will return an ERROR ( findpivot() is called even for the column even though there is no lower column to test if a[n-1][n-1] is zero). If findpivot() returns OK, then uptriangle() calls process_col() to reduce the lower column to zero. We have included debug statements in process_col() to help track the process. The function pr2adbl() prints the two dimensional array - we will soon write this function.

The function findpivot() calls on function findnonzero() to find a non-zero pivot in column k if a[k][k] is zero. If a pivot is found, it swaps the appropriate rows and returns OK. Otherwise, it reurns ERROR. The function findnonzero() merely scans the lower column k for a non-zero element. It either returns the row in which it finds a non-zero element or it returns -1 if no such element is found. Rows of the array are swapped by the function swaprows() which also includes a debug statement to prints the row indices of the rows being swapped.

When uptriangle() returns with OK status, the array will be in upper triangular form. The next step in solving the equations is to employ back substitution to find the values of the unknowns. We now examine the back substitution process. As we saw earlier, we must solve equations backwards starting at index and proceeding to index 0. The equation in upper triangular form looks like this:

a[i][i]*x[i] + a[i][i+1]*x[i+1] + ... + a[i][n-1]*x[n-1] = a[i][n]
Recall, in our representation, the right hand side is the column of the two dimensional array. For each index , we must sum all contributions from those unknowns already solved for, i.e. those x[i] with index greater than . This is the following sum:
sum = a[i][i+1]*x[i+1] + ... + a[i][n-1]*x[n-1]
We then subtract this sum from the right hand side, a[i][n], and divide the result by a[i][i] to determine the solution for x[i]. The algorithm is shown below:
/*   Algorithm: Back_Substitution */
     for (i = n - 1; i >= 0;  i--) {         /* go backwards */
          sum = 0;

for (j = i + 1; j <= n - 1; j++) /* sum all contributions from */ sum += a[i][j] * x[j]; /* x[j] with j > i */ x[i] = (a[i][n] - sum) / a[i][i]; /* solve for x[i] */ }

We can now write the function gauss() that solves a set of equations by the Gaussian elimination method which first calls on uptriangle() to convert the coefficients to upper triangular form. If this succeeds, then back substitution is carried out to find the solutions. As with other functions, gauss() returns OK if successful, and ERROR otherwise. The code is shown in Figure 9.24.

The code is straight forward. It incorporates the back substitution algorithm after the function call to uptriangle(). If the function call returns ERROR, the equations cannot be solved and gauss() returns ERROR. Otherwise, gauss() proceeds with back substitution and stores the result in the array x[]. Since all a[i][i] must be non-zero at this point, we do not really need to test if a[i][i] is zero before using it as a divisor; however, we do so as an added precaution.

We are almost ready to use the function gauss() in a program. Before we can do so; however, we need some utility functions to read and print data. Here are the descriptions of these functions:

getcoeffs():
reads the coefficients and the right hand side values into an array; it returns the number of equations.
pr2adbl():
prints an array with rows and columns.
pr1adbl():
prints a solution array.

All these functions use data of type double. The code is shown in Figure 9.25.

Finally, we are ready to write a program driver as shown in Figure 9.26. The driver first reads coefficients and the right hand side values for a set of equations and then calls on gauss() to solve the equations. During the debug phase, both the original data and the transformed upper triangular version are printed. Finally, if the equations are solved with success, the solution is printed. Otherwise, an error message is printed.

During debugging, the macro DEBUG is defined in gauss.h so that we can track the process. The program loops as long as there are equations to be solved. In each case, it gets coefficients using getcoeffs() and solves them using gauss(). During debug, the program uses pr2adbl() to print the original array and the array after gauss transformation. If the solution is possible, the program prints the solution array using pr1adbl(). Here are several example equation solutions:

Sample Session:

The first two sets of equations are solvable; the last set is not because the second equation in the last set is a multiple of the first equation. Thus these equations are linearly dependent and they cannot be solved uniquely. In this case, after the lower column is reduced to zero, a[1][1] is zero. A pivot is found in row 2, rows 1 and 2 are swapped, and lower column 1 is reduced to zero. However, a[2][2] is now zero, and there is no unique way to solve these equations.

If the coefficients are such that the equations are almost but not quite linearly dependent, the solution can be quite imprecise. An improvement in precision may be obtained by using an element with the largest absolute value as the pivot. Implementation of an improved version of the method is left as an exercise.



Previous: 9.4 Arrays of Pointers
Up: 9.1 Two Dimensional Arrays
Next: 9.6 Common Errors
Previous Page: 9.4 Arrays of Pointers
Next Page: 9.6 Common Errors

tep@wiliki.eng.hawaii.edu
Wed Aug 17 09:20:12 HST 1994