## Iterative Methods

The iterative methods programmed here are designed to be used to solve large sparse systems of linear equations where the direct methods can exceed available machine memory and/or be extremely time-consuming. These methods can also be used to solve systems of linear equations in which the matrices can be stored in available RAM but, in this case, it would be advantageous to modify the routines before using them. The three classical iterative methods are (1) Jacobi's method, (2) Gauss-Seidel's method and (3) the successive overrelaxation method.Given a system of linear equations

**A x**=

**b**, and an initial estimate of the solution,

**x**

^{(0)}which may be taken to be

**0**, Jacobi's method generates a sequence of approximations,

**x**

^{(1)},

**x**

^{(2)}, … by

*x*= [

_{i}^{(k+1)}*b*- Σ

_{i}_{j≠i}

*a*] /

_{ij}x_{j}^{k}*a*, for

_{ii}*i*= 0, … ,

*n*- 1.

**x**

^{(1)},

**x**

^{(2)}, … by

*x*= [

_{i}^{(k+1)}*b*- Σ

_{i}_{j<i}

*a*- Σ

_{ij}x_{j}^{k+1}_{j>i}

*a*] /

_{ij}x_{j}^{k}*a*, for

_{ii}*i*= 0, … ,

*n*- 1.

**x**

^{(1)},

**x**

^{(2)}, … by

*r*= [

_{i}^{(k)}*b*- Σ

_{i}_{j<i}

*a*- Σ

_{ij}x_{j}^{k+1}_{j≥i}

*a*] /

_{ij}x_{j}^{k}*a*, for

_{ii}*i*= 0, … ,

*n*- 1.

*x*=

_{i}^{(k+1)}*x*+ ω

_{i}^{(k)}*r*.

_{i}^{(k)}### Function List

- int Jacobi_Method( void(*iterate)(double*,double*), double x[ ], double x0[ ], int n, double tolerance, int max_tries )

This method is used to solve the linear system of equations**Ax = B**where**A**is an*n×n*matrix,**B**is an*n×1*matrix (or*n*dimensional column vector) for the*n*-dimensional column vector**x**using Jacobi's method. The argument void(*iterate)(double*,double*) is a user supplied function which evaluates**x = [B - (A-D) * x0]D**, where^{-1}**D**is the diagonal of**A**, i.e. (*iterate)(x,x0) uses*x0*and returns*x*via the argument list. The user supplied function should not modify*x0*. The iteration terminates successfully when the max-norm of the difference two successive iterations is within the user specified*tolerance*and terminates unsuccessfully when*max_tries*iterations occur in which the max-norm of the difference of two successive iterations are never within the specified tolerance. The function returns a 0 if it terminates successfully and it returns -1 if the procedure fails to converge within*max_tries*tries.

- int Gauss_Seidel_Method( double (*iterate)(double*,int), double x[ ], double x0[ ], int n, double tolerance, int max_tries )

This method is used to solve the linear system of equations**Ax = B**where**A**is an*n×n*matrix,**B**is an*n×1*matrix (or*n*dimensional column vector) for the*n*-dimensional column vector**x**using the Gauss-Seidel method. The argument double(*iterate)(double*,int) is a user supplied function which evaluates the*i*row of^{th}**x = [B - (A-D) * x0]D**, where^{-1}**D**is the diagonal of**A**, i.e. (*iterate)(x0,i) uses*x0*and returns*x[i]*. Each iteration therefore requires*n*calls to (*iterate). The iteration terminates successfully when the max-norm of the difference two successive iterations is within the user specified*tolerance*and terminates unsuccessfully when*max_tries*iterations occur in which the max-norm of the difference of two successive iterations are never within the specified tolerance. The function returns a 0 if it terminates successfully and it returns -1 if the procedure fails to converge within*max_tries*tries.

- int Successive_Overrelaxation_Method( double (*residual)(double*,int), double x[ ], double x0[ ], int n, double omega, double tolerance, int max_tries )

This method is used to solve the linear system of equations**Ax = B**where**A**is an*n×n*matrix,**B**is an*n×1*matrix (or*n*dimensional column vector) for the*n*-dimensional column vector**x**using the successive overrelaxation method. The argument double(*residual)(x, i) is a user supplied function which evaluates the i-th residual of**{B[i] - A[i][*] x[*]}(1/A[i][i])}**. The next estimate of**x[i]**is the the sum of the old estimate of**x[i]**plus**omega**times the residual returned above. The SOR method converges for**omega**being strictly between 0 and 2. If omega = 1, the method is equivalent to the Gauss-Seidel method.

*C* Source Code

- The file, jacobi.c, contains the routine Jacobi_Method( ).

- The file, gauss_seidel.c, contains the routine Gauss_Seidel_Method( ).

- The file, sor.c, contains the routine Successive_Overrelaxation_Method( ).