Systems of Linear Equations - Lab

Section 1

Untar the package of C code on the class website named 'lineq.tar'. This will create a subdirectory named 'lineq' under the current directory. Change into that 'lineq' subdirectory.

In this section we will use simple Gaussian elimination code to find the solution of two systems of four equations in four unknowns by linking to the functions supplied in the source code file 'lineq.c'. First inspect the source code in 'lineq.c' and its header file 'lineq.h' with

	less lineq.c
	less lineq.h

The first system to solve is:

Prepare to edit the main program C source code file by starting a gedit process running in the background with:

	gedit fourbyfour.c &

Note the ampersand character '&' appended to this command line. This will cause the command to be run in the background but return to a bash shell prompt to allow bash shell commands for compiling and running C code to continue to be run in the terminal. This command starts a gedit process to edit the C source code, and you should see a separate gedit window pop up with the fourbyfour.c source code ready for editing. Read through this short program and see how the main program uses the routine 'alloc_matrix' supplied in the tool file 'lineq.c' to allocate memory for a variable sized two-dimensional matrix, in this case 4 by 4. Find the code near the top of the fourbyfour.c source file and complete the lines indicated with the comments so that the desired matrix corresponding to the system of linear equations, both the coefficient matrix and the right hand side vector, is initialized. Keep in mind these are double precision floating point assignment expressions, so be sure to format your constants appropriately. Also note that although we conventionally write systems of equations with the subscripts on the variables and coefficients starting with 1, in C we must start subscripts to arrays starting with 0!

Save your completed source code file and compile it linked together with the lineq.c code with the shell command:

	gcc -Wall fourbyfour.c lineq.c -o fourbyfour

Note that for this example system we are not calling any functions in the C standard math library, so there is no need to add the '-lm' library search option at the end of this compilation command line (although it would not hurt anything if it were added). Now at the shell prompt, run the compiled program and note the solution it outputs:

	./fourbyfour

Go back to the source code editing window and change the initialization of the right hand side vector so that now the system of equations to be solved is:

Save your edited source code file, recompile and rerun, noting the new solution.

Now we will deliberately give the Gaussian elimination function a system of equations that is indeterminate and has no unique solution. The easiest way to do this is to change one of the equations, say the forth, so that it is simply another of the equations in the system multiplied by a constant. Change the coefficients in the matrix and the right hand side vector entry in the source code to make the last equation be:

Which of the other three equations is this new forth equation a multiple of? Save your edited source code file, recompile and rerun. What lines in the source code are responsible for sensing and reporting this error condition?

Finally see what happens when the memory allocation that sets up space in memory for the coefficient array and vectors representing the system is given an erroneous size. There is one assignment statement in the fourbyfour.c source code that sets the sizes of the square matrix and the vectors that the rest of the program will need. Find that statement, and change the size it specifies from 4 to the erroneous value of 2. Save, recompile and rerun. What happens?

Section 2

Now we will work with a source code file pivot.c that will demonstrate the increased numerical stability that is afforded by adding scaled partial pivoting to the Gaussian elimination algorithm. The equations to be solved are in a simple 2x2 system, where 'd' is some small number:

First verify that this system has a straightforward solution of x₁ = 1/3 and x₂ = 2/3. Prepare to edit the main program C source code file by starting a gedit process running in the background with:

	gedit pivot.c &

Note that in this program we are submitting the same problem to the two solving tools gauss() and gauss_pivotmax() in order to compare their results. Since the first call to gauss() will rewrite entries in the coefficient matrix, we have to reset the coefficient values a second time before the call to gauss_pivotmax(). Also note the additional use of the rindex integer array to hold the row indices that have been arranged by gauss_pivotmax(). These indices must be used at the end to display the solution in the proper order. Find the different spots in the pivot.c source file and complete the lines indicated with the comments so that d is set to a starting value of 1.0e-4 and the desired matrix corresponding to the system of linear equations, both the coefficient matrix and the right hand side vector, is initialized. Keep in mind these are double precision floating point assignment expressions, so be sure to format your constants appropriately. Save your completed source code file and compile it linked together with the lineq.c code with the shell command:

	gcc -Wall pivot.c lineq.c -o pivot

Run your compiled program and note the solution that both algorithms arrive at.

Now edit your program with decreasing values of 'd'. Try 1.0e-6, 1.0e-8, 1.0e-10, 1.0e-12, etc. Recompile and rerun the program for each new value of d. At what value of d do you see numerical instability in the 'naive' gauss() algorithm? How does the gauss_pivotmax() algorithm perform with these same 'd' values?