Gauss–Seidel method

Description

Let $A x = b$ be a square system of n linear equations, where:

$A = [\begin{matrix} a_{11} & a_{12} & \dots & a_{1 n} \\ a_{21} & a_{22} & \dots & a_{2 n} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ a_{n 1} & a_{n 2} & \dots & a_{n n} \end{matrix}], x = [\begin{matrix} x_{1} \\ x_{2} \\ ⋮ \\ x_{n} \end{matrix}], b = [\begin{matrix} b_{1} \\ b_{2} \\ ⋮ \\ b_{n} \end{matrix}] .$

When $A$ and $b$ are known, and $x$ is unknown, the Gauss–Seidel method can be used to iteratively approximate $x$ . The vector $x^{(0)}$ denotes the initial guess for $x$ , often $x_{i}^{(0)} = 0$ for $i = 1, 2, . . ., n$ . Denote by $x^{(k)}$ the $k$ -th approximation or iteration of $x$ , and by $x^{(k + 1)}$ the approximation of $x$ at the next (or $k + 1$ -th) iteration.

Matrix-based formula

The solution is obtained iteratively via $L x^{(k + 1)} = b - U x^{(k)},$ where the matrix $A$ is decomposed into a lower triangular component $L$ , and a strictly upper triangular component $U$ such that $A = L + U$ . More specifically, the decomposition of $A$ into $L$ and $U$ is given by:

$A = \underset{L}{\underset{⏟}{[\begin{matrix} a_{11} & 0 & \dots & 0 \\ a_{21} & a_{22} & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ a_{n 1} & a_{n 2} & \dots & a_{n n} \end{matrix}]}} + \underset{U}{\underset{⏟}{[\begin{matrix} 0 & a_{12} & \dots & a_{1 n} \\ 0 & 0 & \dots & a_{2 n} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & 0 \end{matrix}]}} .$

Why the matrix-based formula works

The system of linear equations may be rewritten as:

\begin{aligned} A x & = b \\ (L + U) x & = b \\ L x + U x & = b \\ L x & = b - U x \end{aligned}

The Gauss–Seidel method now solves the left hand side of this expression for $x$ , using the previous value for $x$ on the right hand side. Analytically, this may be written as $x^{(k + 1)} = L^{- 1} (b - U x^{(k)}) .$

Element-based formula

However, by taking advantage of the triangular form of $L$ , the elements of $x^{(k + 1)}$ can be computed sequentially for each row $i$ using forward substitution: $x_{i}^{(k + 1)} = \frac{1}{a_{i i}} (b_{i} - \sum_{j = 1}^{i - 1} a_{i j} x_{j}^{(k + 1)} - \sum_{j = i + 1}^{n} a_{i j} x_{j}^{(k)}), i = 1, 2, \dots, n .$

Notice that the formula uses two summations per iteration which can be expressed as one summation $\sum_{j \neq i} a_{i j} x_{j}$ that uses the most recently calculated iteration of $x_{j}$ . 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 related to that of the (iterative) Jacobi method, with an important difference:

In Gauss-Seidel, the computation of $x^{(k + 1)}$ uses the elements of $x^{(k + 1)}$ that have already been computed, and only the elements of $x^{(k)}$ that have not been computed in the $(k + 1)$ -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, since they can have a very long critical path, and are thus most feasible for sparse matrices. Furthermore, the values at each iteration are dependent on the order of the original equations.

Gauss-Seidel is the same as successive over-relaxation with $ω = 1$ .

3 more chapters are free with a Lably account.

Convergence

The convergence properties of the Gauss–Seidel method are dependent on the matrix $A$ . Namely, the procedure is known to converge if either:

$A$ is symmetric positive-definite, or
$A$ is strictly or irreducibly diagonally dominant.

The Gauss–Seidel method may converge even if these conditions are not satisfied.

Golub and Van Loan give a theorem for an algorithm that splits $A$ into two parts. Suppose $A = M - N$ is nonsingular. Let $r = ρ (M^{- 1} N)$ be the spectral radius of $M^{- 1} N$ . Then the iterates $x^{(k)}$ defined by $M x^{(k + 1)} = N x^{(k)} + b$ converge to $x = A^{- 1} b$ for any starting vector $x^{(0)}$ if $M$ is nonsingular and $r < 1$ .

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:

<b id="mwlQ">algorithm</b> Gauss–Seidel method <b id="mwlg">is</b> <b id="mwlw">inputs:</b> <var style="padding-right: 1px;" about="#mwt75" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"var","href":"./Template:Var"},"params":{"1":{"wt":"A"}},"i":0}}]}' id="mwmA">A</var>, <var style="padding-right: 1px;" about="#mwt76" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"var","href":"./Template:Var"},"params":{"1":{"wt":"b"}},"i":0}}]}' id="mwmQ">b</var> <span class="nowrap" about="#mwt77" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"nowrap","href":"./Template:Nowrap"},"params":{"1":{"wt":"'''output:''' ''φ''"}},"i":0}}]}' id="mwmg"><b>output:</b> <i>φ</i></span> <span class="nowrap" about="#mwt78" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"nowrap","href":"./Template:Nowrap"},"params":{"1":{"wt":"Choose an initial guess ''φ'' to the solution"}},"i":0}}]}' id="mwmw">Choose an initial guess <i>φ</i> to the solution</span> <b id="mwnA">repeat</b> until convergence <b id="mwnQ">for</b> <var style="padding-right: 1px;" about="#mwt79" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"var","href":"./Template:Var"},"params":{"1":{"wt":"i"}},"i":0}}]}' id="mwng">i</var> <b id="mwnw">from</b> 1 <b id="mwoA">until</b> <var style="padding-right: 1px;" about="#mwt80" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"var","href":"./Template:Var"},"params":{"1":{"wt":"n"}},"i":0}}]}' id="mwoQ">n</var> <b id="mwog">do</b> <span class="nowrap" about="#mwt81" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"nowrap","href":"./Template:Nowrap"},"params":{"1":{"wt":"''σ'' ← 0"}},"i":0}}]}' id="mwow"><i>σ</i> ← 0</span> <b id="mwpA">for</b> <var style="padding-right: 1px;" about="#mwt82" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"var","href":"./Template:Var"},"params":{"1":{"wt":"j"}},"i":0}}]}' id="mwpQ">j</var> <b id="mwpg">from</b> 1 <b id="mwpw">until</b> <var style="padding-right: 1px;" about="#mwt83" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"var","href":"./Template:Var"},"params":{"1":{"wt":"n"}},"i":0}}]}' id="mwqA">n</var> <b id="mwqQ">do</b> <b id="mwqg">if</b> <var style="padding-right: 1px;" about="#mwt84" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"var","href":"./Template:Var"},"params":{"1":{"wt":"j"}},"i":0}}]}' id="mwqw">j</var> ≠ <var style="padding-right: 1px;" about="#mwt85" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"var","href":"./Template:Var"},"params":{"1":{"wt":"i"}},"i":0}}]}' id="mwrA">i</var> <b id="mwrQ">then</b> <span class="nowrap" about="#mwt86" typeof="mw:Transclusion" data-mw="{"parts":[{"template":{"target":{"wt":"nowrap","href":"./Template:Nowrap"},"params":{"1":{"wt":"''σ'' ← ''σ'' + ''a''<sub>''ij''</sub>''φ''<sub>''j''</sub>"}},"i":0}}]}" id="mwrg"><i>σ</i> ← <i>σ</i> + <i>a</i><sub><i>ij</i></sub><i>φ</i><sub><i>j</i></sub></span> <b id="mwrw">end if</b> <b id="mwsA">end</b> (<var style="padding-right: 1px;" about="#mwt87" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"var","href":"./Template:Var"},"params":{"1":{"wt":"j"}},"i":0}}]}' id="mwsQ">j</var>-loop) <span class="nowrap" about="#mwt88" typeof="mw:Transclusion" data-mw="{"parts":[{"template":{"target":{"wt":"nowrap","href":"./Template:Nowrap"},"params":{"1":{"wt":"''φ''<sub>''i''</sub> ← (''b''<sub>''i''</sub> − ''σ'') / ''a''<sub>''ii''</sub>"}},"i":0}}]}" id="mwsg"><i>φ</i><sub><i>i</i></sub> ← (<i>b</i><sub><i>i</i></sub> − <i>σ</i>) / <i>a</i><sub><i>ii</i></sub></span> <b id="mwsw">end</b> (<var style="padding-right: 1px;" about="#mwt89" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"var","href":"./Template:Var"},"params":{"1":{"wt":"i"}},"i":0}}]}' id="mwtA">i</var>-loop) check if convergence is reached <b id="mwtQ">end</b> (repeat)

Examples

An example for the matrix version

A linear system shown as $A x = b$ is given by: $A = [\begin{matrix} 16 & 3 \\ 7 & - 11 \end{matrix}] and b = [\begin{matrix} 11 \\ 13 \end{matrix}] .$

Use the equation $x^{(k + 1)} = L^{- 1} (b - U x^{(k)})$ in the form $x^{(k + 1)} = T x^{(k)} + c$ where:

T = - L^{- 1} U and c = L^{- 1} b .

Decompose $A$ into the sum of a lower triangular component $L$ and a strict upper triangular component $U$ : $L = [\begin{matrix} 16 & 0 \\ 7 & - 11 \end{matrix}] and U = [\begin{matrix} 0 & 3 \\ 0 & 0 \end{matrix}] .$

The inverse of $L$ is: $L^{- 1} = {[\begin{matrix} 16 & 0 \\ 7 & - 11 \end{matrix}]}^{- 1} = [\begin{matrix} 0.0625 & 0.0000 \\ 0.0398 & - 0.0909 \end{matrix}] .$

Now find: $\begin{matrix} T & = - [\begin{matrix} 0.0625 & 0.0000 \\ 0.0398 & - 0.0909 \end{matrix}] [\begin{matrix} 0 & 3 \\ 0 & 0 \end{matrix}] = [\begin{matrix} 0.000 & - 0.1875 \\ 0.000 & - 0.1194 \end{matrix}], \\ c & = [\begin{matrix} 0.0625 & 0.0000 \\ 0.0398 & - 0.0909 \end{matrix}] [\begin{matrix} 11 \\ 13 \end{matrix}] = [\begin{matrix} 0.6875 \\ - 0.7439 \end{matrix}] . \end{matrix}$

With $T$ and $c$ the vectors $x$ can be obtained iteratively.

First of all, choose $x^{(0)}$ , for example $x^{(0)} = [\begin{matrix} 1.0 \\ 1.0 \end{matrix}] .$ The closer the guess to the final solution, the fewer iterations the algorithm will need.

Then calculate: $\begin{matrix} x^{(1)} & = [\begin{matrix} 0.000 & - 0.1875 \\ 0.000 & - 0.1193 \end{matrix}] [\begin{matrix} 1.0 \\ 1.0 \end{matrix}] + [\begin{matrix} 0.6875 \\ - 0.7443 \end{matrix}] = [\begin{matrix} 0.5000 \\ - 0.8636 \end{matrix}] . \\ x^{(2)} & = [\begin{matrix} 0.000 & - 0.1875 \\ 0.000 & - 0.1193 \end{matrix}] [\begin{matrix} 0.5000 \\ - 0.8636 \end{matrix}] + [\begin{matrix} 0.6875 \\ - 0.7443 \end{matrix}] = [\begin{matrix} 0.8494 \\ - 0.6413 \end{matrix}] . \\ x^{(3)} & = [\begin{matrix} 0.000 & - 0.1875 \\ 0.000 & - 0.1193 \end{matrix}] [\begin{matrix} 0.8494 \\ - 0.6413 \end{matrix}] + [\begin{matrix} 0.6875 \\ - 0.7443 \end{matrix}] = [\begin{matrix} 0.8077 \\ - 0.6678 \end{matrix}] . \\ x^{(4)} & = [\begin{matrix} 0.000 & - 0.1875 \\ 0.000 & - 0.1193 \end{matrix}] [\begin{matrix} 0.8077 \\ - 0.6678 \end{matrix}] + [\begin{matrix} 0.6875 \\ - 0.7443 \end{matrix}] = [\begin{matrix} 0.8127 \\ - 0.6646 \end{matrix}] . \\ x^{(5)} & = [\begin{matrix} 0.000 & - 0.1875 \\ 0.000 & - 0.1193 \end{matrix}] [\begin{matrix} 0.8127 \\ - 0.6646 \end{matrix}] + [\begin{matrix} 0.6875 \\ - 0.7443 \end{matrix}] = [\begin{matrix} 0.8121 \\ - 0.6650 \end{matrix}] . \\ x^{(6)} & = [\begin{matrix} 0.000 & - 0.1875 \\ 0.000 & - 0.1193 \end{matrix}] [\begin{matrix} 0.8121 \\ - 0.6650 \end{matrix}] + [\begin{matrix} 0.6875 \\ - 0.7443 \end{matrix}] = [\begin{matrix} 0.8122 \\ - 0.6650 \end{matrix}] . \\ x^{(7)} & = [\begin{matrix} 0.000 & - 0.1875 \\ 0.000 & - 0.1193 \end{matrix}] [\begin{matrix} 0.8122 \\ - 0.6650 \end{matrix}] + [\begin{matrix} 0.6875 \\ - 0.7443 \end{matrix}] = [\begin{matrix} 0.8122 \\ - 0.6650 \end{matrix}] . \end{matrix}$

As expected, the algorithm converges to the solution: $x = A^{- 1} b \approx [\begin{matrix} 0.8122 \\ - 0.6650 \end{matrix}] .$

In fact, the matrix A is strictly diagonally dominant, but not positive definite.

Another example for the matrix version

Another linear system shown as $A x = b$ is given by:

$A = [\begin{matrix} 2 & 3 \\ 5 & 7 \end{matrix}] and b = [\begin{matrix} 11 \\ 13 \end{matrix}] .$

Use the equation $x^{(k + 1)} = L^{- 1} (b - U x^{(k)})$ in the form $x^{(k + 1)} = T x^{(k)} + c$ where:

T = - L^{- 1} U and c = L^{- 1} b .

Decompose $A$ into the sum of a lower triangular component $L$ and a strict upper triangular component $U$ : $L = [\begin{matrix} 2 & 0 \\ 5 & 7 \end{matrix}] and U = [\begin{matrix} 0 & 3 \\ 0 & 0 \end{matrix}] .$

The inverse of $L$ is: $L^{- 1} = {[\begin{matrix} 2 & 0 \\ 5 & 7 \end{matrix}]}^{- 1} = [\begin{matrix} 0.500 & 0.000 \\ - 0.357 & 0.143 \end{matrix}] .$

Now find: $\begin{matrix} T & = - [\begin{matrix} 0.500 & 0.000 \\ - 0.357 & 0.143 \end{matrix}] [\begin{matrix} 0 & 3 \\ 0 & 0 \end{matrix}] = [\begin{matrix} 0.000 & - 1.500 \\ 0.000 & 1.071 \end{matrix}], \\ c & = [\begin{matrix} 0.500 & 0.000 \\ - 0.357 & 0.143 \end{matrix}] [\begin{matrix} 11 \\ 13 \end{matrix}] = [\begin{matrix} 5.500 \\ - 2.071 \end{matrix}] . \end{matrix}$

With $T$ and $c$ the vectors $x$ can be obtained iteratively.

First of all, we have to choose $x^{(0)}$ , for example $x^{(0)} = [\begin{matrix} 1.1 \\ 2.3 \end{matrix}]$

Then calculate: $\begin{matrix} x^{(1)} & = [\begin{matrix} 0 & - 1.500 \\ 0 & 1.071 \end{matrix}] [\begin{matrix} 1.1 \\ 2.3 \end{matrix}] + [\begin{matrix} 5.500 \\ - 2.071 \end{matrix}] = [\begin{matrix} 2.050 \\ 0.393 \end{matrix}] . \\ x^{(2)} & = [\begin{matrix} 0 & - 1.500 \\ 0 & 1.071 \end{matrix}] [\begin{matrix} 2.050 \\ 0.393 \end{matrix}] + [\begin{matrix} 5.500 \\ - 2.071 \end{matrix}] = [\begin{matrix} 4.911 \\ - 1.651 \end{matrix}] . \\ x^{(3)} & = \dots . \end{matrix}$

In a test for convergence we find that the algorithm diverges. In fact, the matrix $A$ is neither diagonally dominant nor positive definite. Then, convergence to the exact solution $x = A^{- 1} b = [\begin{matrix} - 38 \\ 29 \end{matrix}]$ is not guaranteed and, in this case, will not occur.

An example for the equation version

Suppose given $n$ equations and a starting point $x_{0}$ . At any step in a Gauss-Seidel iteration, solve the first equation for $x_{1}$ in terms of $x_{2}, \dots, x_{n}$ ; then solve the second equation for $x_{2}$ in terms of $x_{1}$ just found and the remaining $x_{3}, \dots, x_{n}$ ; and continue to $x_{n}$ . Then, repeat iterations until convergence is achieved, or break if the divergence in the solutions start to diverge beyond a predefined level.

Consider an example: $\begin{array}{rrrrl} 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. \end{array}$

Solving for $x_{1}, x_{2}, x_{3}$ and $x_{4}$ gives: $\begin{matrix} x_{1} & = x_{2} / 10 - x_{3} / 5 + 3 / 5, \\ x_{2} & = x_{1} / 11 + x_{3} / 11 - 3 x_{4} / 11 + 25 / 11, \\ x_{3} & = - x_{1} / 5 + x_{2} / 10 + x_{4} / 10 - 11 / 10, \\ x_{4} & = - 3 x_{2} / 8 + x_{3} / 8 + 15 / 8. \end{matrix}$

Suppose (0, 0, 0, 0) is the initial approximation, then the first approximate solution is given by: $\begin{matrix} 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{matrix}$

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 (1, 2, −1, 1).

An example using Python and NumPy

The following iterative procedure produces the solution vector of a linear system of equations:

<span id="mwATI"></span><span class="kn" id="mwATM">import</span><span class="w" id="mwATQ"> </span><span class="nn" id="mwATU">numpy</span><span class="w" id="mwATY"> </span><span class="k" id="mwATc">as</span><span class="w" id="mwATg"> </span><span class="nn" id="mwATk">np</span> <span class="n" id="mwATo">ITERATION_LIMIT</span> <span class="o" id="mwATs">=</span> <span class="mi" id="mwATw">1000</span> <span class="c1" id="mwAT0"># initialize the matrix</span> <span class="n" id="mwAT4">A</span> <span class="o" id="mwAT8">=</span> <span class="n" id="mwAUA">np</span><span class="o" id="mwAUE">.</span><span class="kp" id="mwAUI">array</span><span class="p" id="mwAUM">(</span> <span class="p" id="mwAUQ">[</span> <span class="p" id="mwAUU">[</span><span class="mf" id="mwAUY">10.0</span><span class="p" id="mwAUc">,</span> <span class="o" id="mwAUg">-</span><span class="mf" id="mwAUk">1.0</span><span class="p" id="mwAUo">,</span> <span class="mf" id="mwAUs">2.0</span><span class="p" id="mwAUw">,</span> <span class="mf" id="mwAU0">0.0</span><span class="p" id="mwAU4">],</span> <span class="p" id="mwAU8">[</span><span class="o" id="mwAVA">-</span><span class="mf" id="mwAVE">1.0</span><span class="p" id="mwAVI">,</span> <span class="mf" id="mwAVM">11.0</span><span class="p" id="mwAVQ">,</span> <span class="o" id="mwAVU">-</span><span class="mf" id="mwAVY">1.0</span><span class="p" id="mwAVc">,</span> <span class="mf" id="mwAVg">3.0</span><span class="p" id="mwAVk">],</span> <span class="p" id="mwAVo">[</span><span class="mf" id="mwAVs">2.0</span><span class="p" id="mwAVw">,</span> <span class="o" id="mwAV0">-</span><span class="mf" id="mwAV4">1.0</span><span class="p" id="mwAV8">,</span> <span class="mf" id="mwAWA">10.0</span><span class="p" id="mwAWE">,</span> <span class="o" id="mwAWI">-</span><span class="mf" id="mwAWM">1.0</span><span class="p" id="mwAWQ">],</span> <span class="p" id="mwAWU">[</span><span class="mf" id="mwAWY">0.0</span><span class="p" id="mwAWc">,</span> <span class="mf" id="mwAWg">3.0</span><span class="p" id="mwAWk">,</span> <span class="o" id="mwAWo">-</span><span class="mf" id="mwAWs">1.0</span><span class="p" id="mwAWw">,</span> <span class="mf" id="mwAW0">8.0</span><span class="p" id="mwAW4">],</span> <span class="p" id="mwAW8">]</span> <span class="p" id="mwAXA">)</span> <span class="c1" id="mwAXE"># initialize the RHS vector</span> <span class="n" id="mwAXI">b</span> <span class="o" id="mwAXM">=</span> <span class="n" id="mwAXQ">np</span><span class="o" id="mwAXU">.</span><span class="kp" id="mwAXY">array</span><span class="p" id="mwAXc">([</span><span class="mf" id="mwAXg">6.0</span><span class="p" id="mwAXk">,</span> <span class="mf" id="mwAXo">25.0</span><span class="p" id="mwAXs">,</span> <span class="o" id="mwAXw">-</span><span class="mf" id="mwAX0">11.0</span><span class="p" id="mwAX4">,</span> <span class="mf" id="mwAX8">15.0</span><span class="p" id="mwAYA">])</span> <span class="nb" id="mwAYE">print</span><span class="p" id="mwAYI">(</span><span class="s2" id="mwAYM">"System of equations:"</span><span class="p" id="mwAYQ">)</span> <span class="k" id="mwAYU">for</span> <span class="n" id="mwAYY">i</span> <span class="ow" id="mwAYc">in</span> <span class="nb" id="mwAYg">range</span><span class="p" id="mwAYk">(</span><span class="n" id="mwAYo">A</span><span class="o" id="mwAYs">.</span><span class="kp" id="mwAYw">shape</span><span class="p" id="mwAY0">[</span><span class="mi" id="mwAY4">0</span><span class="p" id="mwAY8">]):</span> <span class="n" id="mwAZA">row</span> <span class="o" id="mwAZE">=</span> <span class="s2" id="mwAZI">" + "</span><span class="o" id="mwAZM">.</span><span class="n" id="mwAZQ">join</span><span class="p" id="mwAZU">(</span><span class="sa" id="mwAZY">f</span><span class="s2" id="mwAZc">"</span><span class="si" id="mwAZg">{</span><span class="n" id="mwAZk">A</span><span class="p" id="mwAZo">[</span><span class="n" id="mwAZs">i</span><span class="p" id="mwAZw">,</span><span class="w" id="mwAZ0"> </span><span class="n" id="mwAZ4">j</span><span class="p" id="mwAZ8">]</span><span class="si" id="mwAaA">:</span><span class="s2" id="mwAaE">3g</span><span class="si" id="mwAaI">}</span><span class="s2" id="mwAaM">*x</span><span class="si" id="mwAaQ">{</span><span class="n" id="mwAaU">j</span><span class="w" id="mwAaY"> </span><span class="o" id="mwAac">+</span><span class="w" id="mwAag"> </span><span class="mi" id="mwAak">1</span><span class="si" id="mwAao">}</span><span class="s2" id="mwAas">"</span> <span class="k" id="mwAaw">for</span> <span class="n" id="mwAa0">j</span> <span class="ow" id="mwAa4">in</span> <span class="nb" id="mwAa8">range</span><span class="p" id="mwAbA">(</span><span class="n" id="mwAbE">A</span><span class="o" id="mwAbI">.</span><span class="kp" id="mwAbM">shape</span><span class="p" id="mwAbQ">[</span><span class="mi" id="mwAbU">1</span><span class="p" id="mwAbY">]))</span> <span class="nb" id="mwAbc">print</span><span class="p" id="mwAbg">(</span><span class="sa" id="mwAbk">f</span><span class="s2" id="mwAbo">"[</span><span class="si" id="mwAbs">{</span><span class="n" id="mwAbw">row</span><span class="si" id="mwAb0">}</span><span class="s2" id="mwAb4">] = [</span><span class="si" id="mwAb8">{</span><span class="n" id="mwAcA">b</span><span class="p" id="mwAcE">[</span><span class="n" id="mwAcI">i</span><span class="p" id="mwAcM">]</span><span class="si" id="mwAcQ">:</span><span class="s2" id="mwAcU">3g</span><span class="si" id="mwAcY">}</span><span class="s2" id="mwAcc">]"</span><span class="p" id="mwAcg">)</span> <span class="n" id="mwAck">x</span> <span class="o" id="mwAco">=</span> <span class="n" id="mwAcs">np</span><span class="o" id="mwAcw">.</span><span class="kp" id="mwAc0">zeros_like</span><span class="p" id="mwAc4">(</span><span class="n" id="mwAc8">b</span><span class="p" id="mwAdA">)</span> <span class="k" id="mwAdE">for</span> <span class="n" id="mwAdI">it_count</span> <span class="ow" id="mwAdM">in</span> <span class="nb" id="mwAdQ">range</span><span class="p" id="mwAdU">(</span><span class="mi" id="mwAdY">1</span><span class="p" id="mwAdc">,</span> <span class="n" id="mwAdg">ITERATION_LIMIT</span><span class="p" id="mwAdk">):</span> <span class="n" id="mwAdo">x_new</span> <span class="o" id="mwAds">=</span> <span class="n" id="mwAdw">np</span><span class="o" id="mwAd0">.</span><span class="kp" id="mwAd4">zeros_like</span><span class="p" id="mwAd8">(</span><span class="n" id="mwAeA">x</span><span class="p" id="mwAeE">)</span> <span class="nb" id="mwAeI">print</span><span class="p" id="mwAeM">(</span><span class="sa" id="mwAeQ">f</span><span class="s2" id="mwAeU">"Iteration </span><span class="si" id="mwAeY">{</span><span class="n" id="mwAec">it_count</span><span class="si" id="mwAeg">}</span><span class="s2" id="mwAek">: </span><span class="si" id="mwAeo">{</span><span class="n" id="mwAes">x</span><span class="si" id="mwAew">}</span><span class="s2" id="mwAe0">"</span><span class="p" id="mwAe4">)</span> <span class="k" id="mwAe8">for</span> <span class="n" id="mwAfA">i</span> <span class="ow" id="mwAfE">in</span> <span class="nb" id="mwAfI">range</span><span class="p" id="mwAfM">(</span><span class="n" id="mwAfQ">A</span><span class="o" id="mwAfU">.</span><span class="kp" id="mwAfY">shape</span><span class="p" id="mwAfc">[</span><span class="mi" id="mwAfg">0</span><span class="p" id="mwAfk">]):</span> <span class="n" id="mwAfo">s1</span> <span class="o" id="mwAfs">=</span> <span class="n" id="mwAfw">A</span><span class="p" id="mwAf0">[</span><span class="n" id="mwAf4">i</span><span class="p" id="mwAf8">,</span> <span class="p" id="mwAgA">:</span><span class="n" id="mwAgE">i</span><span class="p" id="mwAgI">]</span> <span class="o" id="mwAgM">@</span> <span class="n" id="mwAgQ">x_new</span><span class="p" id="mwAgU">[:</span><span class="n" id="mwAgY">i</span><span class="p" id="mwAgc">]</span> <span class="n" id="mwAgg">s2</span> <span class="o" id="mwAgk">=</span> <span class="n" id="mwAgo">A</span><span class="p" id="mwAgs">[</span><span class="n" id="mwAgw">i</span><span class="p" id="mwAg0">,</span> <span class="n" id="mwAg4">i</span> <span class="o" id="mwAg8">+</span> <span class="mi" id="mwAhA">1</span> <span class="p" id="mwAhE">:]</span> <span class="o" id="mwAhI">@</span> <span class="n" id="mwAhM">x</span><span class="p" id="mwAhQ">[</span><span class="n" id="mwAhU">i</span> <span class="o" id="mwAhY">+</span> <span class="mi" id="mwAhc">1</span> <span class="p" id="mwAhg">:]</span> <span class="n" id="mwAhk">x_new</span><span class="p" id="mwAho">[</span><span class="n" id="mwAhs">i</span><span class="p" id="mwAhw">]</span> <span class="o" id="mwAh0">=</span> <span class="p" id="mwAh4">(</span><span class="n" id="mwAh8">b</span><span class="p" id="mwAiA">[</span><span class="n" id="mwAiE">i</span><span class="p" id="mwAiI">]</span> <span class="o" id="mwAiM">-</span> <span class="n" id="mwAiQ">s1</span> <span class="o" id="mwAiU">-</span> <span class="n" id="mwAiY">s2</span><span class="p" id="mwAic">)</span> <span class="o" id="mwAig">/</span> <span class="n" id="mwAik">A</span><span class="p" id="mwAio">[</span><span class="n" id="mwAis">i</span><span class="p" id="mwAiw">,</span> <span class="n" id="mwAi0">i</span><span class="p" id="mwAi4">]</span> <span class="k" id="mwAi8">if</span> <span class="n" id="mwAjA">np</span><span class="o" id="mwAjE">.</span><span class="kp" id="mwAjI">allclose</span><span class="p" id="mwAjM">(</span><span class="n" id="mwAjQ">x</span><span class="p" id="mwAjU">,</span> <span class="n" id="mwAjY">x_new</span><span class="p" id="mwAjc">,</span> <span class="n" id="mwAjg">rtol</span><span class="o" id="mwAjk">=</span><span class="mf" id="mwAjo">1e-8</span><span class="p" id="mwAjs">):</span> <span class="k" id="mwAjw">break</span> <span class="n" id="mwAj0">x</span> <span class="o" id="mwAj4">=</span> <span class="n" id="mwAj8">x_new</span> <span class="nb" id="mwAkA">print</span><span class="p" id="mwAkE">(</span><span class="sa" id="mwAkI">f</span><span class="s2" id="mwAkM">"Solution: </span><span class="si" id="mwAkQ">{</span><span class="n" id="mwAkU">x</span><span class="si" id="mwAkY">}</span><span class="s2" id="mwAkc">"</span><span class="p" id="mwAkg">)</span> <span class="n" id="mwAkk">error</span> <span class="o" id="mwAko">=</span> <span class="n" id="mwAks">A</span> <span class="o" id="mwAkw">@</span> <span class="n" id="mwAk0">x</span> <span class="o" id="mwAk4">-</span> <span class="n" id="mwAk8">b</span> <span class="nb" id="mwAlA">print</span><span class="p" id="mwAlE">(</span><span class="sa" id="mwAlI">f</span><span class="s2" id="mwAlM">"Error: </span><span class="si" id="mwAlQ">{</span><span class="n" id="mwAlU">error</span><span class="si" id="mwAlY">}</span><span class="s2" id="mwAlc">"</span><span class="p" id="mwAlg">)</span>

Produces the output:

<span id="mwAlw"></span><span class="go" id="mwAl0">System of equations:</span> <span class="go" id="mwAl4">[ 10*x1 + -1*x2 + 2*x3 + 0*x4] = [ 6]</span> <span class="go" id="mwAl8">[ -1*x1 + 11*x2 + -1*x3 + 3*x4] = [ 25]</span> <span class="go" id="mwAmA">[ 2*x1 + -1*x2 + 10*x3 + -1*x4] = [-11]</span> <span class="go" id="mwAmE">[ 0*x1 + 3*x2 + -1*x3 + 8*x4] = [ 15]</span> <span class="go" id="mwAmI">Iteration 1: [ 0. 0. 0. 0.]</span> <span class="go" id="mwAmM">Iteration 2: [ 0.6 2.32727273 -0.98727273 0.87886364]</span> <span class="go" id="mwAmQ">Iteration 3: [ 1.03018182 2.03693802 -1.0144562 0.98434122]</span> <span class="go" id="mwAmU">Iteration 4: [ 1.00658504 2.00355502 -1.00252738 0.99835095]</span> <span class="go" id="mwAmY">Iteration 5: [ 1.00086098 2.00029825 -1.00030728 0.99984975]</span> <span class="go" id="mwAmc">Iteration 6: [ 1.00009128 2.00002134 -1.00003115 0.9999881 ]</span> <span class="go" id="mwAmg">Iteration 7: [ 1.00000836 2.00000117 -1.00000275 0.99999922]</span> <span class="go" id="mwAmk">Iteration 8: [ 1.00000067 2.00000002 -1.00000021 0.99999996]</span> <span class="go" id="mwAmo">Iteration 9: [ 1.00000004 1.99999999 -1.00000001 1. ]</span> <span class="go" id="mwAms">Iteration 10: [ 1. 2. -1. 1.]</span> <span class="go" id="mwAmw">Solution: [ 1. 2. -1. 1.]</span> <span class="go" id="mwAm0">Error: [ 2.06480930e-08 -1.25551054e-08 3.61417563e-11 0.00000000e+00]</span>

Program to solve arbitrary number of equations using Matlab

The following code uses the formula $x_{i}^{(k + 1)} = \frac{1}{a_{i i}} (b_{i} - \sum_{j < i} a_{i j} x_{j}^{(k + 1)} - \sum_{j > i} a_{i j} x_{j}^{(k)}), \begin{array}{l} i = 1, 2, \dots, n \\ k = 0, 1, 2, \dots \end{array}$

<span id="mwAnM"></span><span class="k" id="mwAnQ">function</span><span class="w" id="mwAnU"> </span>x<span class="w" id="mwAnY"> </span><span class="p" id="mwAnc">=</span><span class="w" id="mwAng"> </span><span class="nf" id="mwAnk">gauss_seidel</span><span class="p" id="mwAno">(</span>A, b, x, iters<span class="p" id="mwAns">)</span> <span class="w" id="mwAnw"> </span><span class="k" id="mwAn0">for</span><span class="w" id="mwAn4"> </span><span class="nb" id="mwAn8">i</span><span class="w" id="mwAoA"> </span><span class="p" id="mwAoE">=</span><span class="w" id="mwAoI"> </span><span class="mi" id="mwAoM">1</span><span class="p" id="mwAoQ">:</span><span class="n" id="mwAoU">iters</span> <span class="w" id="mwAoY"> </span><span class="k" id="mwAoc">for</span><span class="w" id="mwAog"> </span><span class="nb" id="mwAok">j</span><span class="w" id="mwAoo"> </span><span class="p" id="mwAos">=</span><span class="w" id="mwAow"> </span><span class="mi" id="mwAo0">1</span><span class="p" id="mwAo4">:</span><span class="nb" id="mwAo8">size</span><span class="p" id="mwApA">(</span><span class="n" id="mwApE">A</span><span class="p" id="mwApI">,</span><span class="mi" id="mwApM">1</span><span class="p" id="mwApQ">)</span> <span class="w" id="mwApU"> </span><span class="n" id="mwApY">x</span><span class="p" id="mwApc">(</span><span class="nb" id="mwApg">j</span><span class="p" id="mwApk">)</span><span class="w" id="mwApo"> </span><span class="p" id="mwAps">=</span><span class="w" id="mwApw"> </span><span class="p" id="mwAp0">(</span><span class="n" id="mwAp4">b</span><span class="p" id="mwAp8">(</span><span class="nb" id="mwAqA">j</span><span class="p" id="mwAqE">)</span><span class="w" id="mwAqI"> </span><span class="o" id="mwAqM">-</span><span class="w" id="mwAqQ"> </span><span class="nb" id="mwAqU">sum</span><span class="p" id="mwAqY">(</span><span class="n" id="mwAqc">A</span><span class="p" id="mwAqg">(</span><span class="nb" id="mwAqk">j</span><span class="p" id="mwAqo">,:)</span><span class="o" id="mwAqs">'.*</span><span class="n" id="mwAqw">x</span><span class="p" id="mwAq0">)</span><span class="w" id="mwAq4"> </span><span class="o" id="mwAq8">+</span><span class="w" id="mwArA"> </span><span class="n" id="mwArE">A</span><span class="p" id="mwArI">(</span><span class="nb" id="mwArM">j</span><span class="p" id="mwArQ">,</span><span class="nb" id="mwArU">j</span><span class="p" id="mwArY">)</span><span class="o" id="mwArc">*</span><span class="n" id="mwArg">x</span><span class="p" id="mwArk">(</span><span class="nb" id="mwAro">j</span><span class="p" id="mwArs">))</span><span class="w" id="mwArw"> </span><span class="o" id="mwAr0">/</span><span class="w" id="mwAr4"> </span><span class="n" id="mwAr8">A</span><span class="p" id="mwAsA">(</span><span class="nb" id="mwAsE">j</span><span class="p" id="mwAsI">,</span><span class="nb" id="mwAsM">j</span><span class="p" id="mwAsQ">);</span> <span class="w" id="mwAsU"> </span><span class="k" id="mwAsY">end</span> <span class="w" id="mwAsc"> </span><span class="k" id="mwAsg">end</span> <span class="k" id="mwAsk">end</span>