Gauss–Seidel method
   HOME

TheInfoList



OR:

In
numerical linear algebra Numerical linear algebra, sometimes called applied linear algebra, is the study of how matrix operations can be used to create computer algorithms which efficiently and accurately provide approximate answers to questions in continuous mathematics. ...
, the Gauss–Seidel method, also known as the Liebmann method or the method of successive displacement, is an
iterative method In computational mathematics, an iterative method is a Algorithm, mathematical procedure that uses an initial value to generate a sequence of improving approximate solutions for a class of problems, in which the ''n''-th approximation is derived fr ...
used to solve a
system of linear equations In mathematics, a system of linear equations (or linear system) is a collection of one or more linear equations involving the same variable (math), variables. For example, :\begin 3x+2y-z=1\\ 2x-2y+4z=-2\\ -x+\fracy-z=0 \end is a system of three ...
. It is named after the
German German(s) may refer to: * Germany (of or related to) **Germania (historical use) * Germans, citizens of Germany, people of German ancestry, or native speakers of the German language ** For citizens of Germany, see also German nationality law **Ger ...
mathematician A mathematician is someone who uses an extensive knowledge of mathematics in their work, typically to solve mathematical problems. Mathematicians are concerned with numbers, data, quantity, structure, space, models, and change. History On ...
s
Carl Friedrich Gauss Johann Carl Friedrich Gauss (; german: Gauß ; la, Carolus Fridericus Gauss; 30 April 177723 February 1855) was a German mathematician and physicist who made significant contributions to many fields in mathematics and science. Sometimes refer ...
and
Philipp Ludwig von Seidel Philipp Ludwig von Seidel (; 24 October 1821 in Zweibrücken, Germany – 13 August 1896 in Munich, German Empire) was a German mathematician. He was the son of Julie Reinhold and Justus Christian Felix Seidel. Lakatos credits von Seidel with di ...
, and is similar to the
Jacobi method In numerical linear algebra, the Jacobi method is an iterative algorithm for determining the solutions of a strictly diagonally dominant system of linear equations. Each diagonal element is solved for, and an approximate value is plugged in. The ...
. Though it can be applied to any matrix with non-zero elements on the diagonals, convergence is only guaranteed if the matrix is either
strictly diagonally dominant In mathematics, a square matrix is said to be diagonally dominant if, for every row of the matrix, the magnitude of the diagonal entry in a row is larger than or equal to the sum of the magnitudes of all the other (non-diagonal) entries in that row ...
, or
symmetric Symmetry (from grc, συμμετρία "agreement in dimensions, due proportion, arrangement") in everyday language refers to a sense of harmonious and beautiful proportion and balance. In mathematics, "symmetry" has a more precise definiti ...
and
positive definite In mathematics, positive definiteness is a property of any object to which a bilinear form or a sesquilinear form may be naturally associated, which is positive-definite. See, in particular: * Positive-definite bilinear form * Positive-definite f ...
. It was only mentioned in a private letter from Gauss to his student
Gerling Gerling were an Australian electronica, alternative rock trio formed in 1993. From early 1997 the members were Darren Cross on guitar and lead vocals, Presser (real name Paul Towner) on drums and Burke Reid on guitar and vocals. Their second al ...
in 1823. A publication was not delivered before 1874 by Seidel.


Description

The Gauss–Seidel method is an iterative technique for solving a square system of linear equations with unknown : A\mathbf x = \mathbf b . It is defined by the iteration L_* \mathbf^ = \mathbf - U \mathbf^, where \mathbf^ is the -th approximation or iteration of \mathbf,\,\mathbf^ is the next or -th iteration of \mathbf, and the matrix is decomposed into a
lower triangular In mathematics, a triangular matrix is a special kind of square matrix. A square matrix is called if all the entries ''above'' the main diagonal are zero. Similarly, a square matrix is called if all the entries ''below'' the main diagonal are ...
component L_*, and a
strictly upper triangular In mathematics, a triangular matrix is a special kind of square matrix. A square matrix is called if all the entries ''above'' the main diagonal are zero. Similarly, a square matrix is called if all the entries ''below'' the main diagonal are ...
component U; that is, A = L_* + U . In more detail, write out , and in their components: A=\begin a_ & a_ & \cdots & a_ \\ a_ & a_ & \cdots & a_ \\ \vdots & \vdots & \ddots & \vdots \\a_ & a_ & \cdots & a_ \end, \qquad \mathbf = \begin x_ \\ x_2 \\ \vdots \\ x_n \end , \qquad \mathbf = \begin b_ \\ b_2 \\ \vdots \\ b_n \end. Then the decomposition of into its lower triangular component and its strictly upper triangular component is given by: A = \underbrace_ + \underbrace_ . The system of linear equations may be rewritten as: :L_* \mathbf = \mathbf - U \mathbf The Gauss–Seidel method now solves the left hand side of this expression for x, using previous value for x on the right hand side. Analytically, this may be written as: \mathbf^ = L_*^ \left(\mathbf - U \mathbf^\right). However, by taking advantage of the triangular form of L_*, the elements of can be computed sequentially using
forward substitution In mathematics, a triangular matrix is a special kind of square matrix. A square matrix is called if all the entries ''above'' the main diagonal are zero. Similarly, a square matrix is called if all the entries ''below'' the main diagonal are ...
: x^_i = \frac \left(b_i - \sum_^a_x^_j - \sum_^a_x^_j \right),\quad i=1,2,\dots,n. The procedure is generally continued until the changes made by an iteration are below some tolerance, such as a sufficiently small residual.


Discussion

The element-wise formula for the Gauss–Seidel method is extremely similar to that of the
Jacobi method In numerical linear algebra, the Jacobi method is an iterative algorithm for determining the solutions of a strictly diagonally dominant system of linear equations. Each diagonal element is solved for, and an approximate value is plugged in. The ...
. The computation of uses the elements of that have already been computed, and only the elements of that have not been computed in the -th iteration. This means that, unlike the Jacobi method, only one storage vector is required as elements can be overwritten as they are computed, which can be advantageous for very large problems. However, unlike the Jacobi method, the computations for each element are generally much harder to implement in
parallel Parallel is a geometric term of location which may refer to: Computing * Parallel algorithm * Parallel computing * Parallel metaheuristic * Parallel (software), a UNIX utility for running programs in parallel * Parallel Sysplex, a cluster of IBM ...
, since they can have a very long critical path, and are thus most feasible for
sparse matrices In numerical analysis and scientific computing, a sparse matrix or sparse array is a matrix in which most of the elements are zero. There is no strict definition regarding the proportion of zero-value elements for a matrix to qualify as sparse b ...
. Furthermore, the values at each iteration are dependent on the order of the original equations. Gauss-Seidel is the same as SOR (successive over-relaxation) with \omega=1.


Convergence

The convergence properties of the Gauss–Seidel method are dependent on the matrix ''A''. Namely, the procedure is known to converge if either: * is symmetric
positive-definite In mathematics, positive definiteness is a property of any object to which a bilinear form or a sesquilinear form may be naturally associated, which is positive-definite. See, in particular: * Positive-definite bilinear form * Positive-definite fu ...
, or * is strictly or irreducibly
diagonally dominant In mathematics, a square matrix is said to be diagonally dominant if, for every row of the matrix, the magnitude of the diagonal entry in a row is larger than or equal to the sum of the magnitudes of all the other (non-diagonal) entries in that row ...
. The Gauss–Seidel method sometimes converges even if these conditions are not satisfied.


Algorithm

Since elements can be overwritten as they are computed in this algorithm, only one storage vector is needed, and vector indexing is omitted. The algorithm goes as follows: Inputs: , repeat until convergence for from 1 until do for from 1 until do if ≠ then end if end (-loop) end (-loop) check if convergence is reached end (repeat)


Examples


An example for the matrix version

A linear system shown as A \mathbf = \mathbf is given by: A= \begin 16 & 3 \\ 7 & -11 \\ \end \quad \text \quad b= \begin 11 \\ 13 \end. We want to use the equation \mathbf^ = L_*^ (\mathbf - U \mathbf^) in the form \mathbf^ = T \mathbf^ + C where: :T = - L_*^ U \quad \text \quad C = L_*^ \mathbf. We must decompose A into the sum of a lower triangular component L_* and a strict upper triangular component U: L_*= \begin 16 & 0 \\ 7 & -11 \\ \end \quad \text \quad U = \begin 0 & 3 \\ 0 & 0 \end. The inverse of L_* is: L_*^ = \begin 16 & 0 \\ 7 & -11 \end^ = \begin 0.0625 & 0.0000 \\ 0.0398 & -0.0909 \\ \end. Now we can find: \begin T &= - \begin 0.0625 & 0.0000 \\ 0.0398 & -0.0909 \end \begin 0 & 3 \\ 0 & 0 \end = \begin 0.000 & -0.1875 \\ 0.000 & -0.1194 \end, \\ ex C &= \begin 0.0625 & 0.0000 \\ 0.0398 & -0.0909 \end \begin 11 \\ 13 \end = \begin 0.6875 \\ -0.7439 \end. \end Now we have T and C and we can use them to obtain the vectors \mathbf iteratively. First of all, we have to choose \mathbf^: we can only guess. The better the guess, the quicker the algorithm will perform. We choose a starting point: x^ = \begin 1.0 \\ 1.0 \end. We can then calculate: \begin x^ &= \begin 0.000 & -0.1875 \\ 0.000 & -0.1193 \end \begin 1.0 \\ 1.0 \end + \begin 0.6875 \\ -0.7443 \end = \begin 0.5000 \\ -0.8636 \end. \\ ex x^ &= \begin 0.000 & -0.1875 \\ 0.000 & -0.1193 \end \begin 0.5000 \\ -0.8636 \end + \begin 0.6875 \\ -0.7443 \end = \begin 0.8494 \\ -0.6413 \end. \\ ex x^ &= \begin 0.000 & -0.1875 \\ 0.000 & -0.1193 \end \begin 0.8494 \\ -0.6413 \\ \end + \begin 0.6875 \\ -0.7443 \end = \begin 0.8077 \\ -0.6678 \end. \\ ex x^ &= \begin 0.000 & -0.1875 \\ 0.000 & -0.1193 \end \begin 0.8077 \\ -0.6678 \end + \begin 0.6875 \\ -0.7443 \end = \begin 0.8127 \\ -0.6646 \end. \\ ex x^ &= \begin 0.000 & -0.1875 \\ 0.000 & -0.1193 \end \begin 0.8127 \\ -0.6646 \end + \begin 0.6875 \\ -0.7443 \end = \begin 0.8121 \\ -0.6650 \end. \\ ex x^ &= \begin 0.000 & -0.1875 \\ 0.000 & -0.1193 \end \begin 0.8121 \\ -0.6650 \end + \begin 0.6875 \\ -0.7443 \end = \begin 0.8122 \\ -0.6650 \end. \\ ex x^ &= \begin 0.000 & -0.1875 \\ 0.000 & -0.1193 \end \begin 0.8122 \\ -0.6650 \end + \begin 0.6875 \\ -0.7443 \end = \begin 0.8122 \\ -0.6650 \end. \end As expected, the algorithm converges to the exact solution: \mathbf = A^ \mathbf \approx \begin 0.8122\\ -0.6650 \end. In fact, the matrix is strictly diagonally dominant (but not positive definite).


Another example for the matrix version

Another linear system shown as A \mathbf = \mathbf is given by: A= \begin 2 & 3 \\ 5 & 7 \\ \end \quad \text \quad b = \begin 11 \\ 13 \\ \end. We want to use the equation \mathbf^ = L_*^ (\mathbf - U \mathbf^) in the form \mathbf^ = T \mathbf^ + C where: :T = - L_*^ U \quad \text \quad C = L_*^ \mathbf. We must decompose A into the sum of a lower triangular component L_* and a strict upper triangular component U: L_*= \begin 2 & 0 \\ 5 & 7 \\ \end \quad \text \quad U = \begin 0 & 3 \\ 0 & 0 \\ \end. The inverse of L_* is: L_*^ = \begin 2 & 0 \\ 5 & 7 \\ \end^ = \begin 0.500 & 0.000 \\ -0.357 & 0.143 \\ \end . Now we can find: \begin T &= - \begin 0.500 & 0.000 \\ -0.357 & 0.143 \\ \end \begin 0 & 3 \\ 0 & 0 \\ \end = \begin 0.000 & -1.500 \\ 0.000 & 1.071 \\ \end, \\ ex C &= \begin 0.500 & 0.000 \\ -0.357 & 0.143 \\ \end \begin 11 \\ 13 \\ \end = \begin 5.500 \\ -2.071 \\ \end. \end Now we have T and C and we can use them to obtain the vectors \mathbf iteratively. First of all, we have to choose \mathbf^: we can only guess. The better the guess, the quicker will perform the algorithm. We suppose: x^ = \begin 1.1 \\ 2.3 \end. We can then calculate: \begin x^ &= \begin 0 & -1.500 \\ 0 & 1.071 \\ \end \begin 1.1 \\ 2.3 \\ \end + \begin 5.500 \\ -2.071 \\ \end = \begin 2.050 \\ 0.393 \\ \end. \\ ex x^ &= \begin 0 & -1.500 \\ 0 & 1.071 \\ \end \begin 2.050 \\ 0.393 \\ \end + \begin 5.500 \\ -2.071 \\ \end = \begin 4.911 \\ -1.651 \end. \\ ex x^ &= \cdots. \end If we test for convergence we'll find that the algorithm diverges. In fact, the matrix A is neither diagonally dominant nor positive definite. Then, convergence to the exact solution \mathbf = A^ \mathbf = \begin -38\\ 29 \end is not guaranteed and, in this case, will not occur.


An example for the equation version

Suppose given equations where ''x''''n'' are vectors of these equations and starting point ''x''0. From the first equation solve for ''x''1 in terms of x_, x_, \dots, x_n. For the next equations substitute the previous values of ''x''s. To make it clear consider an example. \begin 10x_1 &- x_2 &+ 2x_3 & & = 6, \\ -x_1 &+ 11x_2 &- x_3 &+ 3x_4 & = 25, \\ 2x_1 &- x_2 &+ 10x_3 &- x_4 & = -11, \\ & 3x_2 &- x_3 &+ 8x_4 & = 15. \end Solving for x_1, x_2, x_3 and x_4 gives: \begin x_1 & = x_2/10 - x_3/5 + 3/5, \\ x_2 & = x_1/11 + x_3/11 - 3x_4/11 + 25/11, \\ x_3 & = -x_1/5 + x_2/10 + x_4/10 - 11/10, \\ x_4 & = -3x_2/8 + x_3/8 + 15/8. \end Suppose we choose as the initial approximation, then the first approximate solution is given by \begin x_1 & = 3/5 = 0.6, \\ x_2 & = (3/5)/11 + 25/11 = 3/55 + 25/11 = 2.3272, \\ x_3 & = -(3/5)/5 +(2.3272)/10-11/10 = -3/25 + 0.23272-1.1 = -0.9873,\\ x_4 & = -3(2.3272)/8 +(-0.9873)/8+15/8 = 0.8789. \end Using the approximations obtained, the iterative procedure is repeated until the desired accuracy has been reached. The following are the approximated solutions after four iterations. The exact solution of the system is .


An example using Python and NumPy

The following numerical procedure simply iterates to produce the solution vector. import numpy as np ITERATION_LIMIT = 1000 # initialize the matrix A = np.array(
10., -1., 2., 0. 1 (one, unit, unity) is a number representing a single or the only entity. 1 is also a numerical digit and represents a single unit (measurement), unit of counting or measurement. For example, a line segment of ''unit length'' is a line segment ...
1., 11., -1., 3. ., -1., 10., -1. ., 3., -1., 8.) # initialize the RHS vector b = np.array( .0, 25.0, -11.0, 15.0 print("System of equations:") for i in range(A.shape : row = "*x"_for_j_in_range(A.shape[1.html" ;"title=".html" ;"title=""*x" for j in range(A.shape[1">"*x" for j in range(A.shape[1">.html" ;"title=""*x" for j in range(A.shape[1">"*x" for j in range(A.shape[1 print("[] = []".format(" + ".join(row), b[i])) x = np.zeros_like(b) for it_count in range(1, ITERATION_LIMIT): x_new = np.zeros_like(x) print(f"Iteration : ") for i in range(A.shape : s1 = np.dot(A , :i x_new i s2 = np.dot(A , i + 1 : x + 1 : x_new = (b - s1 - s2) / A
, i The comma is a punctuation mark that appears in several variants in different languages. It has the same shape as an apostrophe or single closing quotation mark () in many typefaces, but it differs from them in being placed on the baseline ...
if np.allclose(x, x_new, rtol=1e-8): break x = x_new print(f"Solution: ") error = np.dot(A, x) - b print(f"Error: ")
Produces the output: System of equations: 10*x1 + -1*x2 + 2*x3 + 0*x4= 6 -1*x1 + 11*x2 + -1*x3 + 3*x4= 25 2*x1 + -1*x2 + 10*x3 + -1*x4= 11 0*x1 + 3*x2 + -1*x3 + 8*x4= 15Iteration 1: 0. 0. 0. 0.Iteration 2: 0.6 2.32727273 -0.98727273 0.87886364Iteration 3:
1.03018182 2.03693802 -1.0144562 0.98434122 1. is the ordinal form of the number one in a number of European languages. 1. may also refer to: Association football Austria * 1. FC Vöcklabruck, a defunct Austrian association football club Czech Republic * 1. FC Karlovy Vary, a Czech asso ...
Iteration 4: 1.00658504 2.00355502 -1.00252738 0.99835095Iteration 5:
1.00086098 2.00029825 -1.00030728 0.99984975 1. is the ordinal form of the number one in a number of European languages. 1. may also refer to: Association football Austria * 1. FC Vöcklabruck, a defunct Austrian association football club Czech Republic * 1. FC Karlovy Vary, a Czech asso ...
Iteration 6: 1.00009128 2.00002134 -1.00003115 0.9999881 Iteration 7: 1.00000836 2.00000117 -1.00000275 0.99999922Iteration 8: 1.00000067 2.00000002 -1.00000021 0.99999996Iteration 9:
1.00000004 1.99999999 -1.00000001 1. 1. is the ordinal form of the number one in a number of European languages. 1. may also refer to: Association football Austria * 1. FC Vöcklabruck, a defunct Austrian association football club Czech Republic * 1. FC Karlovy Vary, a Czech asso ...
Iteration 10:
1. 2. -1. 1. 1. is the Ordinal number (linguistics), ordinal form of the number 1 (number), one in a number of Languages of Europe, European languages. 1. may also refer to: Association football Austria * 1. FC Vöcklabruck, a defunct Austrian association foot ...
Solution:
1. 2. -1. 1. 1. is the Ordinal number (linguistics), ordinal form of the number 1 (number), one in a number of Languages of Europe, European languages. 1. may also refer to: Association football Austria * 1. FC Vöcklabruck, a defunct Austrian association foot ...
Error: 2.06480930e-08 -1.25551054e-08 3.61417563e-11 0.00000000e+00


Program to solve arbitrary no. of equations using Matlab

The following code uses the formula x^_i = \frac \left(b_i - \sum_a_x^_j - \sum_a_x^_j \right),\quad \begin i=1,2,\ldots,n \\ k=0,1,2,\ldots \end function x = gauss_seidel(A, b, x, iters) for i = 1:iters for j = 1:size(A,1) x(j) = (b(j) - sum(A(j,:)'.*x) + A(j,j)*x(j)) / A(j,j); end end end


See also

* Gaussian belief propagation * Iterative method. Linear systems *
Kaczmarz method The Kaczmarz method or Kaczmarz's algorithm is an iterative algorithm for solving linear equation systems A x = b . It was first discovered by the Polish mathematician Stefan Kaczmarz, and was rediscovered in the field of image reconstruction from ...
(a "row-oriented" method, whereas Gauss-Seidel is "column-oriented." See, for example
this paper
) *
Matrix splitting In the mathematical discipline of numerical linear algebra, a matrix splitting is an expression which represents a given matrix as a sum or difference of matrices. Many iterative methods (for example, for systems of differential equations) depen ...
*
Richardson iteration Modified Richardson iteration is an iterative method for solving a system of linear equations. Richardson iteration was proposed by Lewis Fry Richardson in his work dated 1910. It is similar to the Jacobi and Gauss–Seidel method. We seek the ...


Notes


References

* . * . *


External links

*
Gauss–Seidel from www.math-linux.com
From Holistic Numerical Methods Institute
BicksonMatlab code
{{DEFAULTSORT:Gauss-Seidel Method Numerical linear algebra Articles with example pseudocode Relaxation (iterative methods) Articles with example Python (programming language) code Articles with example MATLAB/Octave code