#define MAX_DIM 100 #define MAX_ITER 500 #define TOLERANCE 1.e-6 void gauss_seid

Question

#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 error = 1.0;
    double ite_error;
    double e, sum, temp;
    int i, j;

    while (error > 1.e-6)
    {
        ite_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 > ite_error)
            ite_error = e;
    }
    error = ite_error;
}

Related Questions

Navigate

• Browse (All)
• Browse #
• Subjects

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.