Thank you for visiting our site! You landed on this page because you entered a search term similar to this: Solving Non Linear Simultaneous Equations, here's the result:
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, occurin a variety of mathematical applications in science, engineering, economics,and social sciences. Examples include: electronic circuit analysis, econometricanalysis, structural analysis, etc. In the most general case, the number ofequations, , may be different from the number of unknowns,
;thus, it may not bepossible to find a unique solution. However, if
equals
, there is a goodchance of finding a unique solution for the unknowns.
Our next task is to solve a set of linear algebraic equations, assumingthat the number of equations equals the number of unknowns:
LINEQNS: Read the coefficients and the right hand side values for a set oflinear equations; solve the equations for the unknowns.
The solution of a set of linear equations is fairly complex. We will firstreview the process of solution and then develop the algorithm in small parts.As we develop parts of the algorithm, we will implement these parts asfunctions. The driver will just read the coefficients, call on a function tosolve the equations, and call a function to print the solution.
Let us start with an example of a set of three simultaneous equations in threeunknowns: ,
, and
.
We can use arraysto represent such a set of equations; a two dimensional arrayto store the coefficients, a one dimensional array to store the valuesof the unknowns when solved, and another one dimensional array to store thevalues on the right hand side. Later, we will include the right hand sidevalues as an additional column in the two dimensional array of coefficients.Each row of the two dimensional array stores the coefficients of one of theequations. Since the array index in C starts at 0, we will assume the unknownsare the elements x[0], x[1], and x[2].Similarly, the elements in row zero arethe coefficients in the equation number 0, the elements in row one are forequation 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 onedimensional arrays: x[0], x[1],
A common method for solving such equations is called the Gaussian elimination method. Themethod eliminates (i.e. makes zero) all coefficients below the main diagonal ofthe two dimensional array. It does so by adding multiples of some equations toothers in a systematic way. The elimination makes the array of new coefficientshave an upper triangular form since the lower triangular coefficients are allzero.
The modified equivalent set of equations in
unknowns in the uppertriangular 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. Backsubstitution first solves the last equation which has only oneunknown, x[n-1].It is easily solved for this value --- x[n-1] = y[n-1]/a[n-1][n-1].The next tothe last equation may then be solved -since x[n-1] has been determined already, this value is substituted inthe equation, andthis equation has again only one unknown, x[n-2].The unknown, x[n-2], is solved for,and the process continuesbackward to the next higher equation. At each stage, the values of the unknownssolved for in the previous equations are substituted in the new equationleaving only one unknown. In this manner, each equation has only one unknownwhich 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] = 3We first reduce to zero the coefficients in the firstcolumn below the main diagonal (i.e. array index zero).If the first equation is multiplied by -2 andadded to the second equation, the coefficient in the second row and firstcolumn 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] = 3Similarly, if the first equation is multiplied by -1 and added to the thirdequation, 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] = -3Coefficients 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 multipliedby a multiplier and added to equations below the second; thus, multiplying thesecond equation by -2 and adding to the third makes the coefficient in thesecond 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] = 9We now have equivalent equations with an upper triangular form for the non-zerocoefficients. 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 andsolving for x[1] gives us x[1] = 1.Finally, substituting x[2] and x[1] in thefirst equation gives us x[0] = 1.
From the above discussion, we can see thata general algorithm involves two steps: modify the coefficients of theequations to an upper triangular form, and solve the equations by backsubstitution.
Let us first consider the process of modifyingthe equations to an upper triangularform. Since only the coefficients and the right hand side values are involvedin the computations that modify the equations to upper triangular form, we canwork 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 columnsbelow the main diagonal, storing the modified new values of the elementsin the same elements of the array.Now, it is time to reduce the
lower columnto zero (by lower column, we mean the part of the column below the maindiagonal). 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 weneed to do is exchange this
row with any higher indexedrow 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 triangularform: for each row find a pivot, and reduce the corresponding lower column tozero. If a non-zero pivot element is not found, then one or more equations arelinear combinations of others, the equations are called linearly dependent, and they cannot be solved.
Figures 9.22 and 9.23 showthe set of functions that convert the first rows and columns of anarray to an upper triangular form. These and other functions use a user definedtype, 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, sothe two dimensional array must have MAXrows and MAX+1 columns.Figure 9.21 includes the header file with thedefines and function prototypes used in the program.
Since precision is important in these computations, we have used formalparameters of type double.The two dimensional arrays can store coefficientsfor a maximum of MAX equations (rows) and have MAX + 1columns to accommodatethe right hand side values.
The function uptriangle() transforms coefficientsof the equations to an upper triangular form.For each k from 0 through n-1, itcalls 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 eventhough there is no lower
column to testif 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 andreturns OK.Otherwise, it reurns ERROR.The function findnonzero() merely scans the lowercolumn k for a non-zero element.It either returns the row in which it finds anon-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 statementto 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 substitutionto find the values of the unknowns.We now examine the back substitution process.As we saw earlier, we must solve equationsbackwards starting at index and proceeding to index 0.The
equation inupper 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 thensubtract this sum from the right hand side, a[i][n],and divide the resultby 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 solvesa set of equations by the Gaussian eliminationmethod which first calls on uptriangle() to convertthe coefficients to uppertriangular form. If this succeeds, then back substitution is carried out to findthe solutions.As with other functions, gauss() returns OKif successful, and ERROR otherwise.The code is shown in Figure 9.24.
The code is straight forward.It incorporates the back substitutionalgorithm 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 backsubstitution and stores the result in the array x[].Since all a[i][i] must benon-zero at this point, we do not really need to test if a[i][i]is zero beforeusing 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 doso; however, we need some utility functions to read and print data.Here are thedescriptions of these functions: getcoeffs():
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 coefficientsand 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 thetransformed upper triangular version are printed. Finally, if the equations aresolved with success, the solution is printed. Otherwise, an error message isprinted.
During debugging, the macro DEBUG is defined in gauss.h so that we can track the process.The programloops as long as there are equations to be solved. In each case, it getscoefficients using getcoeffs() and solves them using gauss(). Duringdebug, the program uses pr2adbl() to print the originalarray and the arrayafter gauss transformation. If the solution is possible, the program prints thesolution 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 thesecond equation in the last set is a multiple of the first equation.Thus theseequations are linearly dependent and they cannot be solved uniquely. In thiscase, after the lower column is reduced to zero, a[1][1]is zero. Apivot is found in row 2, rows 1 and 2 are swapped, and lower column 1 isreduced 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 quitelinearly dependent, the solution can be quite imprecise. An improvement inprecision may be obtained by using an element with the largest absolute value as thepivot. Implementation of an improved version of the method is left as anexercise.
Previous Page: 9.4 Arrays of Pointers
Next Page: 9.6 Common Errors
Wed Aug 17 09:20:12 HST 1994