Menu
Home
People
Places
Arts
History
Plants & Animals
Science
Life & Culture
Technology
Reference.org
Jacobi method
open-in-new
<h2 id="description">Description</h2> <p>Let A x = b {\displaystyle A\mathbf {x} =\mathbf {b} } be a square system of <i>n</i> linear equations, where: A = [ a 11 a 12 ⋯ a 1 n a 21 a 22 ⋯ a 2 n ⋮ ⋮ ⋱ ⋮ a n 1 a n 2 ⋯ a n n ] , x = [ x 1 x 2 ⋮ x n ] , b = [ b 1 b 2 ⋮ b n ] . {\displaystyle A={\begin{bmatrix}a_{11}&a_{12}&\cdots &a_{1n}\\a_{21}&a_{22}&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &a_{nn}\end{bmatrix}},\qquad \mathbf {x} ={\begin{bmatrix}x_{1}\\x_{2}\\\vdots \\x_{n}\end{bmatrix}},\qquad \mathbf {b} ={\begin{bmatrix}b_{1}\\b_{2}\\\vdots \\b_{n}\end{bmatrix}}.} </p><p>When A {\displaystyle A} and b {\displaystyle \mathbf {b} } are known, and x {\displaystyle \mathbf {x} } is unknown, we can use the Jacobi method to approximate x {\displaystyle \mathbf {x} } . The vector x ( 0 ) {\displaystyle \mathbf {x} ^{(0)}} denotes our initial guess for x {\displaystyle \mathbf {x} } (often x i ( 0 ) = 0 {\displaystyle \mathbf {x} _{i}^{(0)}=0} for i = 1 , 2 , . . . , n {\displaystyle i=1,2,...,n} ). We denote x ( k ) {\displaystyle \mathbf {x} ^{(k)}} as the <i>k-</i>th approximation or iteration of x {\displaystyle \mathbf {x} } , and x ( k + 1 ) {\displaystyle \mathbf {x} ^{(k+1)}} is the next (or <i>k</i>+1) iteration of x {\displaystyle \mathbf {x} } . </p> <h3>Matrix-based formula</h3> <p>Then <i>A</i> can be decomposed into a <a href="/facts/Diagonal_matrix/tjIAHrE6">diagonal</a> component <i>D</i>, a lower triangular part <i>L</i> and an upper triangular part <i>U</i>: A = D + L + U where D = [ a 11 0 ⋯ 0 0 a 22 ⋯ 0 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ a n n ] and L + U = [ 0 a 12 ⋯ a 1 n a 21 0 ⋯ a 2 n ⋮ ⋮ ⋱ ⋮ a n 1 a n 2 ⋯ 0 ] . {\displaystyle A=D+L+U\qquad {\text{where}}\qquad D={\begin{bmatrix}a_{11}&0&\cdots &0\\0&a_{22}&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &a_{nn}\end{bmatrix}}{\text{ and }}L+U={\begin{bmatrix}0&a_{12}&\cdots &a_{1n}\\a_{21}&0&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &0\end{bmatrix}}.} The solution is then obtained iteratively via </p> x ( k + 1 ) = D − 1 ( b − ( L + U ) x ( k ) ) . {\displaystyle \mathbf {x} ^{(k+1)}=D^{-1}(\mathbf {b} -(L+U)\mathbf {x} ^{(k)}).} <h3>Element-based formula</h3> <p>The element-based formula for each row i {\displaystyle i} is thus: x i ( k + 1 ) = 1 a i i ( b i − ∑ j ≠ i a i j x j ( k ) ) , i = 1 , 2 , … , n . {\displaystyle x_{i}^{(k+1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j\neq i}a_{ij}x_{j}^{(k)}\right),\quad i=1,2,\ldots ,n.} The computation of x i ( k + 1 ) {\displaystyle x_{i}^{(k+1)}} requires each element in x ( k ) {\displaystyle \mathbf {x} ^{(k)}} except itself. Unlike the <a href="/facts/Gauss%25E2%2580%2593Seidel_method/eK94QeFt">Gauss–Seidel method</a>, we cannot overwrite x i ( k ) {\displaystyle x_{i}^{(k)}} with x i ( k + 1 ) {\displaystyle x_{i}^{(k+1)}} , as that value will be needed by the rest of the computation. The minimum amount of storage is two vectors of size <i>n</i>. </p> <h2 id="algorithm">Algorithm</h2> Input: initial guess <i>x</i>(0) to the solution, (diagonal dominant) matrix <i>A</i>, right-hand side vector <i>b</i>, convergence criterion Output: solution when convergence is reached Comments: pseudocode based on the element-based formula above <i>k</i> = 0 while convergence not reached do for <i>i</i> := 1 step until n do <i>σ</i> = 0 for <i>j</i> := 1 step until n do if <i>j</i> ≠ <i>i</i> then <i>σ</i> = <i>σ</i> + <i>a</i><i>ij</i> <i>x</i><i>j</i>(<i>k</i>) end end <i>x</i><i>i</i>(<i>k</i>+1) = (<i>b</i><i>i</i> − <i>σ</i>) / <i>a</i><i>ii</i> end increment <i>k</i> end <h2 id="convergence">Convergence</h2> <p>The standard convergence condition (for any iterative method) is when the <a href="/facts/Spectral_radius/j8T2xrLu">spectral radius</a> of the iteration matrix is less than 1: </p> ρ ( D − 1 ( L + U ) ) < 1. {\displaystyle \rho (D^{-1}(L+U))<1.} <p>A sufficient (but not necessary) condition for the method to converge is that the matrix <i>A</i> is strictly or irreducibly <a href="/facts/Diagonally_dominant_matrix/qbRQV7CV">diagonally dominant</a>. Strict row diagonal dominance means that for each row, the absolute value of the diagonal term is greater than the sum of absolute values of other terms: </p> | a i i | > ∑ j ≠ i | a i j | . {\displaystyle \left|a_{ii}\right|>\sum _{j\neq i}{\left|a_{ij}\right|}.} <p>The Jacobi method sometimes converges even if these conditions are not satisfied. </p><p>Note that the Jacobi method does not converge for every symmetric <a href="/facts/Positive-definite_matrix/Hbr6AuPS">positive-definite matrix</a>. For example, A = ( 29 2 1 2 6 1 1 1 1 5 ) ⇒ D − 1 ( L + U ) = ( 0 2 29 1 29 1 3 0 1 6 5 5 0 ) ⇒ ρ ( D − 1 ( L + U ) ) ≈ 1.0661 . {\displaystyle A={\begin{pmatrix}29&2&1\\2&6&1\\1&1&{\frac {1}{5}}\end{pmatrix}}\quad \Rightarrow \quad D^{-1}(L+U)={\begin{pmatrix}0&{\frac {2}{29}}&{\frac {1}{29}}\\{\frac {1}{3}}&0&{\frac {1}{6}}\\5&5&0\end{pmatrix}}\quad \Rightarrow \quad \rho (D^{-1}(L+U))\approx 1.0661\,.} </p> <h2 id="examples">Examples</h2> <h3>Example question</h3> <p>A linear system of the form A x = b {\displaystyle Ax=b} with initial estimate x ( 0 ) {\displaystyle x^{(0)}} is given by </p> A = [ 2 1 5 7 ] , b = [ 11 13 ] and x ( 0 ) = [ 1 1 ] . {\displaystyle A={\begin{bmatrix}2&1\\5&7\\\end{bmatrix}},\ b={\begin{bmatrix}11\\13\\\end{bmatrix}}\quad {\text{and}}\quad x^{(0)}={\begin{bmatrix}1\\1\\\end{bmatrix}}.} <p>We use the equation x ( k + 1 ) = D − 1 ( b − ( L + U ) x ( k ) ) {\displaystyle x^{(k+1)}=D^{-1}(b-(L+U)x^{(k)})} , described above, to estimate x {\displaystyle x} . First, we rewrite the equation in a more convenient form D − 1 ( b − ( L + U ) x ( k ) ) = T x ( k ) + C {\displaystyle D^{-1}(b-(L+U)x^{(k)})=Tx^{(k)}+C} , where T = − D − 1 ( L + U ) {\displaystyle T=-D^{-1}(L+U)} and C = D − 1 b {\displaystyle C=D^{-1}b} . From the known values D − 1 = [ 1 / 2 0 0 1 / 7 ] , L = [ 0 0 5 0 ] and U = [ 0 1 0 0 ] . {\displaystyle D^{-1}={\begin{bmatrix}1/2&0\\0&1/7\\\end{bmatrix}},\ L={\begin{bmatrix}0&0\\5&0\\\end{bmatrix}}\quad {\text{and}}\quad U={\begin{bmatrix}0&1\\0&0\\\end{bmatrix}}.} we determine T = − D − 1 ( L + U ) {\displaystyle T=-D^{-1}(L+U)} as T = [ 1 / 2 0 0 1 / 7 ] { [ 0 0 − 5 0 ] + [ 0 − 1 0 0 ] } = [ 0 − 1 / 2 − 5 / 7 0 ] . {\displaystyle T={\begin{bmatrix}1/2&0\\0&1/7\\\end{bmatrix}}\left\{{\begin{bmatrix}0&0\\-5&0\\\end{bmatrix}}+{\begin{bmatrix}0&-1\\0&0\\\end{bmatrix}}\right\}={\begin{bmatrix}0&-1/2\\-5/7&0\\\end{bmatrix}}.} Further, C {\displaystyle C} is found as C = [ 1 / 2 0 0 1 / 7 ] [ 11 13 ] = [ 11 / 2 13 / 7 ] . {\displaystyle C={\begin{bmatrix}1/2&0\\0&1/7\\\end{bmatrix}}{\begin{bmatrix}11\\13\\\end{bmatrix}}={\begin{bmatrix}11/2\\13/7\\\end{bmatrix}}.} With T {\displaystyle T} and C {\displaystyle C} calculated, we estimate x {\displaystyle x} as x ( 1 ) = T x ( 0 ) + C {\displaystyle x^{(1)}=Tx^{(0)}+C} : x ( 1 ) = [ 0 − 1 / 2 − 5 / 7 0 ] [ 1 1 ] + [ 11 / 2 13 / 7 ] = [ 5.0 8 / 7 ] ≈ [ 5 1.143 ] . {\displaystyle x^{(1)}={\begin{bmatrix}0&-1/2\\-5/7&0\\\end{bmatrix}}{\begin{bmatrix}1\\1\\\end{bmatrix}}+{\begin{bmatrix}11/2\\13/7\\\end{bmatrix}}={\begin{bmatrix}5.0\\8/7\\\end{bmatrix}}\approx {\begin{bmatrix}5\\1.143\\\end{bmatrix}}.} The next iteration yields x ( 2 ) = [ 0 − 1 / 2 − 5 / 7 0 ] [ 5.0 8 / 7 ] + [ 11 / 2 13 / 7 ] = [ 69 / 14 − 12 / 7 ] ≈ [ 4.929 − 1.714 ] . {\displaystyle x^{(2)}={\begin{bmatrix}0&-1/2\\-5/7&0\\\end{bmatrix}}{\begin{bmatrix}5.0\\8/7\\\end{bmatrix}}+{\begin{bmatrix}11/2\\13/7\\\end{bmatrix}}={\begin{bmatrix}69/14\\-12/7\\\end{bmatrix}}\approx {\begin{bmatrix}4.929\\-1.714\\\end{bmatrix}}.} This process is repeated until convergence (i.e., until ‖ A x ( n ) − b ‖ {\displaystyle \|Ax^{(n)}-b\|} is small). The solution after 25 iterations is </p> x = [ 7.111 − 3.222 ] . {\displaystyle x={\begin{bmatrix}7.111\\-3.222\end{bmatrix}}.} <h3>Example question 2</h3> <p>Suppose we are given the following linear system: </p> 10 x 1 − x 2 + 2 x 3 = 6 , − x 1 + 11 x 2 − x 3 + 3 x 4 = 25 , 2 x 1 − x 2 + 10 x 3 − x 4 = − 11 , 3 x 2 − x 3 + 8 x 4 = 15. {\displaystyle {\begin{aligned}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{aligned}}} <p>If we choose (0, 0, 0, 0) as the initial approximation, then the first approximate solution is given by x 1 = ( 6 + 0 − ( 2 ∗ 0 ) ) / 10 = 0.6 , x 2 = ( 25 + 0 + 0 − ( 3 ∗ 0 ) ) / 11 = 25 / 11 = 2.2727 , x 3 = ( − 11 − ( 2 ∗ 0 ) + 0 + 0 ) / 10 = − 1.1 , x 4 = ( 15 − ( 3 ∗ 0 ) + 0 ) / 8 = 1.875. {\displaystyle {\begin{aligned}x_{1}&=(6+0-(2*0))/10=0.6,\\x_{2}&=(25+0+0-(3*0))/11=25/11=2.2727,\\x_{3}&=(-11-(2*0)+0+0)/10=-1.1,\\x_{4}&=(15-(3*0)+0)/8=1.875.\end{aligned}}} Using the approximations obtained, the iterative procedure is repeated until the desired accuracy has been reached. The following are the approximated solutions after five iterations. </p> <table><tbody><tr><th> x 1 {\displaystyle x_{1}} </th><th> x 2 {\displaystyle x_{2}} </th><th> x 3 {\displaystyle x_{3}} </th><th> x 4 {\displaystyle x_{4}} </th></tr><tr><td>0.6</td><td>2.27272</td><td>-1.1</td><td>1.875</td></tr><tr><td>1.04727</td><td>1.7159</td><td>-0.80522</td><td>0.88522</td></tr><tr><td>0.93263</td><td>2.05330</td><td>-1.0493</td><td>1.13088</td></tr><tr><td>1.01519</td><td>1.95369</td><td>-0.9681</td><td>0.97384</td></tr><tr><td>0.98899</td><td>2.0114</td><td>-1.0102</td><td>1.02135</td></tr></tbody></table> <p>The exact solution of the system is (1, 2, −1, 1). </p> <h3>Python example</h3> import numpy as np ITERATION_LIMIT = 1000 # initialize the matrix A = np.array([[10., -1., 2., 0.], [-1., 11., -1., 3.], [2., -1., 10., -1.], [0.0, 3., -1., 8.]]) # initialize the RHS vector b = np.array([6., 25., -11., 15.]) # prints the system print("System:") for i in range(A.shape[0]): row = [f"{A[i, j]}*x{j + 1}" for j in range(A.shape[1])] print(f'{" + ".join(row)} = {b[i]}') print() x = np.zeros_like(b) for it_count in range(ITERATION_LIMIT): if it_count != 0: print(f"Iteration {it_count}: {x}") x_new = np.zeros_like(x) for i in range(A.shape[0]): s1 = np.dot(A[i, :i], x[:i]) s2 = np.dot(A[i, i + 1:], x[i + 1:]) x_new[i] = (b[i] - s1 - s2) / A[i, i] if x_new[i] == x_new[i-1]: break if np.allclose(x, x_new, atol=1e-10, rtol=0.): break x = x_new print("Solution: ") print(x) error = np.dot(A, x) - b print("Error:") print(error) <h2 id="weighted-jacobi-method">Weighted Jacobi method</h2> <p>The weighted Jacobi iteration uses a parameter ω {\displaystyle \omega } to compute the iteration as </p> x ( k + 1 ) = ω D − 1 ( b − ( L + U ) x ( k ) ) + ( 1 − ω ) x ( k ) {\displaystyle \mathbf {x} ^{(k+1)}=\omega D^{-1}(\mathbf {b} -(L+U)\mathbf {x} ^{(k)})+\left(1-\omega \right)\mathbf {x} ^{(k)}} <p>with ω = 2 / 3 {\displaystyle \omega =2/3} being the usual choice.<a class="footnote-ref" id="fnref:1" href="#fn:1"><sup>1</sup></a> From the relation L + U = A − D {\displaystyle L+U=A-D} , this may also be expressed as </p> x ( k + 1 ) = ω D − 1 b + ( I − ω D − 1 A ) x ( k ) {\displaystyle \mathbf {x} ^{(k+1)}=\omega D^{-1}\mathbf {b} +\left(I-\omega D^{-1}A\right)\mathbf {x} ^{(k)}} . <h3>Convergence in the symmetric positive definite case</h3> <p>In case that the system matrix A {\displaystyle A} is of symmetric <a href="/facts/Positive-definite_matrix/Hbr6AuPS">positive-definite</a> type one can show convergence. </p><p>Let C = C ω = I − ω D − 1 A {\displaystyle C=C_{\omega }=I-\omega D^{-1}A} be the iteration matrix. Then, convergence is guaranteed for </p> ρ ( C ω ) < 1 ⟺ 0 < ω < 2 λ max ( D − 1 A ) , {\displaystyle \rho (C_{\omega })<1\quad \Longleftrightarrow \quad 0<\omega <{\frac {2}{\lambda _{\text{max}}(D^{-1}A)}}\,,} <p>where λ max {\displaystyle \lambda _{\text{max}}} is the maximal eigenvalue. </p><p>The spectral radius can be minimized for a particular choice of ω = ω opt {\displaystyle \omega =\omega _{\text{opt}}} as follows min ω ρ ( C ω ) = ρ ( C ω opt ) = 1 − 2 κ ( D − 1 A ) + 1 for ω opt := 2 λ min ( D − 1 A ) + λ max ( D − 1 A ) , {\displaystyle \min _{\omega }\rho (C_{\omega })=\rho (C_{\omega _{\text{opt}}})=1-{\frac {2}{\kappa (D^{-1}A)+1}}\quad {\text{for}}\quad \omega _{\text{opt}}:={\frac {2}{\lambda _{\text{min}}(D^{-1}A)+\lambda _{\text{max}}(D^{-1}A)}}\,,} where κ {\displaystyle \kappa } is the <a href="/facts/Condition_number/9Tjtk6wh">matrix condition number</a>. </p> <h2 id="see-also">See also</h2> <ul><li><a href="/facts/Gauss%25E2%2580%2593Seidel_method/eK94QeFt">Gauss–Seidel method</a></li> <li><a href="/facts/Successive_over-relaxation/rCUV11er">Successive over-relaxation</a></li> <li><a href="/facts/Iterative_method/uKMAJhEh">Iterative method § Linear systems</a></li> <li><a href="/facts/Belief_propagation/TRY5LVbP">Gaussian Belief Propagation</a></li> <li><a href="/facts/Matrix_splitting/lSExxJBj">Matrix splitting</a></li></ul> <h2 id="external-links">External links</h2> <ul><li>This article incorporates text from the article <a href="http://www.cfd-online.com/Wiki/Jacobi_method">Jacobi_method</a> on <a href="http://www.cfd-online.com/Wiki/Main_Page">CFD-Wiki</a> that is under the <a href="/facts/GFDL/IhVHF3zP">GFDL</a> license.</li></ul> <ul><li>Black, Noel; Moore, Shirley & Weisstein, Eric W. <a href="https://mathworld.wolfram.com/JacobiMethod.html">"Jacobi method"</a>. <i><a href="/facts/MathWorld/yc1jodbq">MathWorld</a></i>.</li> <li><a href="http://www.math-linux.com/spip.php?article49">Jacobi Method from www.math-linux.com</a></li></ul> <h2 id="references">References</h2> <ol> <li id="fn:1"><p>Saad, Yousef (2003). Iterative Methods for Sparse Linear Systems (2nd ed.). SIAM. p. 414. ISBN 0898715342. <a href="0898715342" target="_blank">0898715342</a> <a href="#fnref:1" class="footnote-back-ref">↩</a></p></li> </ol>