Part 1: Complete the program given below that implements the Gauss-Seidel iterat
ID: 3683628 • Letter: P
Question
Part 1:
Complete the program given below that implements the Gauss-Seidel iterative method for a general system of n equations with n unknowns ( n should be an input variable). Limit the maximum number of equations that can be solved at any one time to 100. Your code should be able to solve the matrix equation Ax = b , and print out the vector x . The program:
? Reads in the size of the system to be solved (n).
? Then, it reads in all the elements of the system matrix A and of the right-hand-side vector b .
? Next, it tests whether the sufficient convergence criterion for the method is met.
? Next, it solves the system by calling a function that implements the Gauss-Seidel method.
You need to complete the function. Recall that the method is iterative and will not
always converge; therefore, you should include a maximum number of allowed iterations.
? Test your code with the following system:
2x1 - x2 = 1
-x1 + 3x2 - 2x3 = 1.5
-2x2 + 5x3 - 3x4 = 2.5
-3x3 + 3x4 = 1.5
Use an error tolerance of 1.e-6, and the following initial estimates x1(0) = x2(0) = x3(0) = x4(0) = 0. Recall that the Gauss-Seidel method is based on the following solution updates:
where the superscript k denotes the k -th update of an unknown xi , Aij are the system’s coefficients (the elements of the system matrix A ), and n is the dimension of the system (total number of unknowns). In the above expression, conventional indexing, typically used in mathematical sciences, is employed for A, x, and b, i.e., from 1 to n (as opposed to C-indexing that goes from 0 to n-1).
Complete the Program below:
/* code to solve a nxn system using the Gauss-Seidel method */
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define MAX_DIM 100
#define MAX_ITER 500
#define TOLERANCE 1.e-6
void gauss_seidel (double a[][MAX_DIM], double b[], double x[], int n);
void main ()
{
int i, j, n;
int violation_counter, answer;
int violation_rows[MAX_DIM];
double sum;
double a[MAX_DIM][MAX_DIM];
double b[MAX_DIM],x[MAX_DIM];
/* read in data */
n = MAX_DIM + 1;
while (n > MAX_DIM) {
printf ("Enter the dimension of the system to be solved:");
scanf ("%d",&n);
}
printf(" ");
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
printf ("Enter a[%d][%d] element of the system matrix: ",i,j);
scanf ("%lf",&a[i][j]);
}
printf ("Enter b[%d] of the right-hand-side vector: ",i);
scanf ("%lf",&b[i]);
printf (" ");
}
/* test the convergence criterion */
violation_counter = 0;
for (i = 0; i < n; i++) {
sum = 0.0;
for (j = 0; j < n; j++)
if (i != j)
sum = sum + fabs(a[i][j]);
if (fabs(a[i][i]) < sum) {
violation_rows[violation_counter]= i;
violation_counter = violation_counter + 1;
}
if (a[i][i] == 0.0) {
printf ("Found diagonal element equal to zero;
rearrange equations; exiting ... ");
exit (0);
}
}
if (violation_counter > 0) {
printf ("The Gauss-Seidel convergence criterion is violated
in %d rows out of %d ", violation_counter, n);
printf ("Specifically, it was violated in rows: ");
for (i = 0; i < violation_counter; i++)
printf("%d ",violation_rows[i]);
printf(" ");
printf ("Enter 1 if you want to continue; any other
number to abort: ");
scanf ("%d",&answer);
if (answer != 1)
exit(1);
printf ("Check results carefully ");
}
/* initialize the solution vector -- initial guesses */
for (i = 0; i < n; i++) {
printf ("Enter an initial guess for x[%d] of the
solution vector: ",i);
scanf ("%lf",&x[i]);
}
/* solve the system */
gauss_seidel (a,b,x,n);
/* output solution */
for (i = 0; i < n; i++)
printf ("x[%d]=%f ",i,x[i]);
printf(" ");
}
/* function to solve a system using Gauss-Seidel */
void gauss_seidel (double a[][MAX_DIM], double b[], double x[], int n)
{
Explanation / Answer
void gauss_seidel (double a[][MAX_DIM], double b[], double x[], int n)
{
double maxerror = 1.0;
double iteration_error;
double e, sum, temp;
int i, j;
while (maxerror > 1.e-6)
{
iteration_error = 0.0;
for (i = 0; i < n; i++)
{
sum = 0.0;
for (j = 0; j < n; j++)
{
if (i != j)
sum = sum + (a[i][j] * x[j]);
}
}
temp = (a[i][n] - sum)/a[i][i];
e = fabs((temp - x[i]) / x[i]);
x[i]= temp;
if (e > iteration_error)
iteration_error = e;
}
maxerror = iteration_error;
}