# Systems of Linear Equations

 Home Matrix Home Systems of Linear Equations Upper / Lower Triagular Matrices Gaussian Elimination Equilibrate Matrix Choleski Decomposition Crout Decomposition Doolittle Decomposition Tridiagonal Systems Iterative Methods Singular Value Decomposition

## 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
xi(k+1) = [bi - Σj≠i aij xjk] / aii, for i = 0, … , n - 1.
Gauss-Seidel's method generates a sequence of approximations, x(1), x(2), … by
xi(k+1) = [bi - Σj<i aij xjk+1 - Σj>i aij xjk] / aii, for i = 0, … , n - 1.
And the successive overrelaxation method generates a sequence of approximations, x(1), x(2), … by
ri(k) = [bi - Σj<i aij xjk+1 - Σj≥i aij xjk] / aii, for i = 0, … , n - 1.
xi(k+1) = xi(k) + ω ri(k).
The parameter ω is called the relaxation parameter.

### 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-1, where 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 ith row of x = [B - (A-D) * x0]D-1, where 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( ).