Fourteenth Week Examples -- Iterative Methods    Solving n×n Systems of Equations    Finding Largest Eigenvalue

Jacobi and Gauss-Seidel Iteration

```/***********************************************************************************
Treibergs                                                              Mar. 24, 2006

Jacobi Iteration and Gauss-Seidel Iteration Methods

For large matrices it sometimes more efficient to use an iterative scheme to solve a
linear system  Ax = b.

We assume that the matrix is diagonal dominant.

Then  replace A = D + N  where  D  is the diagonal matrix and  N = A - D. Then the
equation may be rewritten

Dx = b - Nx

Then the Jacobi method is to iterate this equation. Thus the vectors
x_1, x_2, x_3, ... tend to the solution

Dx_(k+1) = b - N x_k.

The Gauss-Seidel method is to use the updated  x  as soon as it comes available.
e.g. if the size = 3,

x_(k+1)[1]  =
(b[1] - a[1][2] * x_k[2] - a[1][3] * x_k[3] ) / a[1][1]

x_(k+1)[2]  =
(b[2] - a[2][1] * x_(k+1)[1] - a[2][3] * x_k [3] ) / a[1][1]

x_(k+1)[3]  =
(b[3] - a[3][1] * x_(k+1)[1] - a[3][2] * x_(k+1)[2] ) / a[1][1]

Gauss-Seidel converges faster than the Jacobi method.

A measure of the rate of convergence is the condition number, which measures the
diagonal dominance of the matrix.  For the Jacobi method the condition number is

M = max_(1 <= i <= n) {  sum_{j=1,...,i-1,i+1,...,n} |a[i][j]| / |a[i][i]|  }

If M<1, it guarantees that the Jacobi Method converges at an exponential rate
| soln - x_k | <= M^k | soln - x_1|.
The method may converge for other condition numbers, and there are sharper
convergence tests. The Gauss-Seidel method sometimes converges even if the Jacobi
method does not.

jagase.c
*************************************************************************************/

# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# define  MAXDIM 20
# define TOL 0.5e-7

typedef struct
{
double m[MAXDIM][MAXDIM];
int size;
} matrix;

typedef struct
{
double m[MAXDIM];
int size;
} vector;

matrix
enterm ( void );

vector
enterv ( matrix a );

vector
jacite ( matrix a, vector b, vector v );

vector
gausei ( matrix a, vector b, vector v );

void
printm ( matrix , vector );

void
printv( vector );

int
main ( void )
{
double m=0.0, sum;
int i,j;
matrix a;
vector b,v, w;

a = enterm ();
b = enterv ( a );
printf("\n\nInitial Matrix and Right Side");
printm ( a, b );
for (i = 0; i<b.size; i++)
{
sum=0.0;
v.m[i] = b.m[i]/a.m[i][i];
for(j=0; j< a.size; j++)
if ( j != i)
sum += fabs(a.m[i][j]);
sum /= fabs(a.m[i][i]);
if( sum > m)
m = sum;
}
printf("The condition number of the matrix is  %f\n",m);
v.size=b.size;
w=v;

printf("\n\nJacobi Iterates\n");
printv( v );
for ( i = 0; i < 20; i++)
{
v = jacite( a, b, v);
printv ( v );
}

printf("\n\nGauss-Seidel Iterates\n");
printv( w );
for ( i = 0; i < 20; i++)
{
w = gausei( a, b, w);
printv ( w );
}

return EXIT_SUCCESS;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

matrix
enterm ( void )
{
int i=0, j, ir, ic;
matrix g;

while ( i++ <= 10 )
{
printf ( "\nEnter the size = number of rows = number of columns  : " );
scanf ( "%d", &j );

if ( (1<=j) && (j <= MAXDIM)  )
{
for ( ir=0; ir< j; ir++ )
{
printf ( "Enter row  a(%d,*) :", ir+1 );
for ( ic = 0; ic < j; ic++ )
scanf ( "%lf", &g.m[ir][ic] );
}
g.size = j;
return g;
}
else
printf ( "Bad input: size must be between %d and %d.\n", 1, MAXDIM );

}
exit ( 0 );
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

void
printm ( matrix a, vector v )
{
int i, j;
printf ( "\n\n" );
for ( i = 0; i < a.size; i++ )
{
for ( j = 0; j < a.size; j++ )
printf ( "\t%f", a.m[i][j] );
printf ( "\t  %f\n\n",v.m[i] );
}
return;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

vector
enterv ( matrix a )
{
int i;
vector v;

printf ( "\n\nEnter right side  (b[*])^T  : " );
for ( i = 0; i < a.size; i++ )
scanf ( "%lf", &v.m[i] );
v.size=a.size;
printf("\n");
return v;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

void
printv ( vector v )
{
int j;

for ( j = 0; j < v.size; j++ )
printf ( "\t%f", v.m[j] );
printf ( "\n" );

return;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

vector
jacite ( matrix a, vector b, vector v )
{
double s;
int i, j;
vector c;

for ( i = 0; i < a.size; i++ )
{
s = b.m[i];

for ( j = 0; j < a.size; j++ )
{
if( i != j)
s -=  a.m[i][j] * v.m[j];

}
c.m[i]=s/a.m[i][i];
}
c.size=a.size;
return c;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

vector
gausei ( matrix a, vector b, vector v )
{
double s;
int i, j;
vector c;

for ( i = 0; i < a.size; i++ )
{
s = b.m[i];

for ( j = 0; j < a.size; j++ )
{
if( j < i)
s -=  a.m[i][j] * c.m[j];
else if (j>i)
s -=  a.m[i][j] * v.m[j];
}
c.m[i]=s/a.m[i][i];
}
c.size=a.size;
return c;
}

```

Power Method to find the Largest Eigenvalue

```/***********************************************************************************
Treibergs                                                              Mar. 24, 2006

Power method to find a the largest eigenvalue.

An eigenvalue of a square  n x n  martrix  A  is a complex number  c  for which there
is a nonzero vector  v  such that  Av = cv.  A number is an eigenvalue if and only
if the number satisfies the  n-th degree polynomial equation

p(c)  =  det(A - cI)  =  0,

where  I  is the identity matrix.  p(c)  is called the characteristic polynomial.
It has degree  n  so that there are at most  n  distinct eigenvalues for any given
matrix.  One way to find  eigenvalues is to solve for  c  in  p(c) = 0,  but this
is difficult. Assuming that the eigenvalue of largest magnitude is real, which is
often true of matrices that come up in practice, such as real symmetric matrices,
and that if there is only one eigenvalue with this largest magnitue, which is
true for most matrices, then there is a very simple way to find the eigenvalue
using the power method. The basic idea is to consider the sequence

x_0  any nonzero starting vector.  Then for all k=1,2,3... we just multiply by  A
repeatedly

w_{k+1} x_{k+i} =   A x_k

where w_k is a scaling factor, say, so that x_k does not blow up. For example
letting  w_k  be the length of  A x_k  works, in effext, x_{k+1}  is confined to be
a unit length vector. Under these hypotheses,  w_k  tends to the largest
eigenvalue of  A.

The program asks for the matrix, It iterates  n  times. It prints the vector  x_k
and the eigenvalue  w_k  each loop.

power.c
*************************************************************************************/

# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# define  MAXDIM 20
# define TOL 0.5e-7

typedef struct
{
double m[MAXDIM][MAXDIM];
int size;
} matrix;

typedef struct
{
double m[MAXDIM];
int size;
} vector;

matrix
enterm ( void );

vector
mmultv ( matrix a,  vector v );

vector
scalv ( double c, vector w );

double
vdotv ( vector v, vector w );

void
printm ( matrix  );

void
printv( vector );

int
main ( void )
{
double c;
int i, j, n = 25;
matrix a;
vector v;

a = enterm ();
printf ( "\n\nInitial Matrix" );
printm ( a );

v.size = a.size;
for (i = 0; i < a.size; i++ )         /* pick a starting vextor x_0  */
v.m[i] = fabs(a.m[i][i])+0.0001;

printf ( "\n\nPower Method\n"
"   n        Eigenvalue         Eigenvector\n                           ");
c = sqrt(vdotv( v,v));
v = scalv ( 1.0/c, v );
printv ( v );

for ( i = 1; i <= n; i++)
{                               /* Apply  A  and divide result by its length  */
v = mmultv( a, v);
c = sqrt(vdotv( v,v));
printf ( "\n%4d %20.15f", i, c );
v = scalv ( 1.0/c, v );
printv ( v );
}

printf ( "\n" );
return EXIT_SUCCESS;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

matrix
enterm ( void )
{
int i=0, j, ir, ic;
matrix g;

while ( i++ <= 10 )
{
printf ( "\nEnter the size = number of rows = number of columns  : " );
scanf ( "%d", &j );

if ( (1<=j) && (j <= MAXDIM)  )
{
for ( ir=0; ir< j; ir++ )
{
printf ( "Enter row  a(%d,*) :", ir+1 );
for ( ic = 0; ic < j; ic++ )
scanf ( "%lf", &g.m[ir][ic] );
}
g.size = j;
return g;
}
else
printf ( "Bad input: size must be between %d and %d.\n", 1, MAXDIM );

}
exit ( 0 );
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

void
printm ( matrix a)
{
int i, j;
printf ( "\n\n" );
for ( i = 0; i < a.size; i++ )
{
for ( j = 0; j < a.size; j++ )
printf ( "\t%f", a.m[i][j] );
printf ( "\n\n" );
}
return;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

void
printv ( vector v )
{
int j;
for ( j = 0; j < v.size; j++ )
printf ( "\t%f", v.m[j] );
return;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

vector
mmultv ( matrix a,  vector v )
{
int i,k;
double sum;
vector w;
for(i = 0; i < a.size; i++)
{
sum = 0.0;
for(k = 0; k < a.size; k++)
sum +=  a.m[i][k]*v.m[k];
w.m[i]=sum;
}
w.size=a.size;
return w;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

vector
scalv ( double c, vector v )
{
int k;
vector w;
for(k = 0; k < v.size; k++)
w.m[k]=  c*v.m[k];

w.size=v.size;
return w;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

double
vdotv ( vector v, vector w )
{
int i,j,k;
double sum;
sum = 0.0;
for(k = 0; k < v.size; k++)
sum +=  w.m[k]*v.m[k];
return sum;
}

```