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
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
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],
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 -a[i][k] / a[k][k]The
/* 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
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:
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.