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;
}