Successive over-relaxation

It was devised simultaneously by David M. Young Jr. and by Stanley P. Frankel in 1950 for the purpose of automatically solving linear systems on digital computers. Over-relaxation methods had been used before the work of Young and Frankel. An example is the method of Lewis Fry Richardson, and the methods developed by R. V. Southwell. However, these methods were designed for computation by human calculators, requiring some expertise to ensure convergence to the solution which made them inapplicable for programming on digital computers. These aspects are discussed in the thesis of David M. Young Jr.

Formulation

Given a square system of n linear equations with unknown x:

$A x = b$

where:

$\begin{matrix} 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}] . \end{matrix}$

Then A can be decomposed into a diagonal component D, and strictly lower and upper triangular components L and U:

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

The system of linear equations may be rewritten as:

$(D + ω L) x = ω b - [ω U + (ω - 1) D] x$

for a constant ω > 1, called the relaxation factor.

The method of successive over-relaxation is an iterative technique that 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:

$\begin{matrix} x^{(k + 1)} & = {(D + ω L)}^{- 1} (ω b - [ω U + (ω - 1) D] x^{(k)}) \\ = L_{ω} x^{(k)} + c, \end{matrix}$

where $x^{(k)}$ is the kth approximation or iteration of $x$ and $x^{(k + 1)}$ is the next or k + 1 iteration of $x$ . However, by taking advantage of the triangular form of (D+ωL), the elements of x^(k+1) can be computed sequentially using forward substitution:

$x_{i}^{(k + 1)} = (1 - ω) x_{i}^{(k)} + \frac{ω}{a_{i i}} (b_{i} - \sum_{j < i} a_{i j} x_{j}^{(k + 1)} - \sum_{j > i} a_{i j} x_{j}^{(k)}), i = 1, 2, \dots, n .$

This can again be written analytically in matrix-vector form without the need of inverting the matrix $(D + ω L)$ :

$x^{(k + 1)} = (1 - ω) x^{(k)} + ω D^{- 1} (b - L x^{(k + 1)} - U x^{(k)}) .$

5 more chapters are free with a Lably account.

Convergence

The choice of relaxation factor ω is not necessarily easy, and depends upon the properties of the coefficient matrix. In 1947, Ostrowski proved that if $A$ is symmetric and positive-definite then $ρ (L_{ω}) < 1$ for $0 < ω < 2$ . Thus, convergence of the iteration process follows, but we are generally interested in faster convergence rather than just convergence.

Convergence Rate

The convergence rate for the SOR method can be analytically derived. One needs to assume the following

the relaxation parameter is appropriate: $ω \in (0, 2)$
Jacobi's iteration matrix $C_{Jac} := I - D^{- 1} A$ has only real eigenvalues
Jacobi's method is convergent: $μ := ρ (C_{Jac}) < 1$
the matrix decomposition $A = D + L + U$ satisfies the property that $det (λ D + z L + \frac{1}{z} U) = det (λ D + L + U)$ for any $z \in C ∖ {0}$ and $λ \in C$ .

Then the convergence rate can be expressed as $ρ (C_{ω}) = {\begin{cases} \frac{1}{4} {[ω μ + \sqrt{ω^{2} μ^{2} - 4 (ω - 1)}]}^{2}, & 0 < ω \leq ω_{opt} \\ ω - 1, & ω_{opt} < ω < 2 \end{cases}$ where the optimal relaxation parameter is given by $ω_{opt} := 1 + {(\frac{μ}{1 + \sqrt{1 - μ^{2}}})}^{2} = 1 + \frac{μ^{2}}{4} + O (μ^{3}) .$ In particular, for $ω = 1$ (Gauss-Seidel) it holds that $ρ (C_{ω}) = μ^{2} = ρ (C_{Jac})^{2}$ . For the optimal $ω$ we get $ρ (C_{ω}) = \frac{1 - \sqrt{1 - μ^{2}}}{1 + \sqrt{1 - μ^{2}}} = \frac{μ^{2}}{4} + O (μ^{3})$ , which shows SOR is roughly four times more efficient than Gauss–Seidel.

The last assumption is satisfied for tridiagonal matrices since $Z (λ D + L + U) Z^{- 1} = λ D + z L + \frac{1}{z} U$ for diagonal $Z$ with entries $Z_{i i} = z^{i - 1}$ and $det (λ D + L + U) = det (Z (λ D + L + U) Z^{- 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:

Inputs: <span class="texhtml mvar" style="font-style:italic;" about="#mwt50" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"mvar","href":"./Template:Mvar"},"params":{"1":{"wt":"A"}},"i":0}}]}' id="mwiA">A</span>, <span class="texhtml mvar" style="font-style:italic;" about="#mwt51" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"mvar","href":"./Template:Mvar"},"params":{"1":{"wt":"b"}},"i":0}}]}' id="mwiQ">b</span>, <span class="texhtml mvar" style="font-style:italic;" about="#mwt52" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"mvar","href":"./Template:Mvar"},"params":{"1":{"wt":"ω"}},"i":0}}]}' id="mwig">ω</span> Output: <span class="mwe-math-element mwe-math-element-inline" typeof="mw:Extension/math" about="#mwt53" data-mw='{"name":"math","attrs":{},"body":{"extsrc":"\\phi"}}' id="mwiw"><span class="mwe-math-mathml-a11y" style="display: none;"><math xmlns="http://www.w3.org/1998/Math/MathML" class="mathjax_ignore" alttext="{\displaystyle \phi }"> <semantics> <mrow class="MJX-TeXAtom-ORD"> <mstyle displaystyle="true" scriptlevel="0"> <mi>ϕ<!-- ϕ --></mi> </mstyle> </mrow> <annotation encoding="application/x-tex">{\displaystyle \phi }</annotation> </semantics> </math></span><img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/72b1f30316670aee6270a28334bdf4f5072cdde4" class="mwe-math-fallback-image-inline mw-invert skin-invert" aria-hidden="true" style="vertical-align: -0.671ex; width:1.385ex; height:2.509ex;" alt="{\displaystyle \phi }"/></span> <!-- the leading whitespace in the next line is significant --> Choose an initial guess <span class="mwe-math-element mwe-math-element-inline" typeof="mw:Extension/math" about="#mwt54" data-mw='{"name":"math","attrs":{},"body":{"extsrc":"\\phi"}}' id="mwjA"><span class="mwe-math-mathml-a11y" style="display: none;"><math xmlns="http://www.w3.org/1998/Math/MathML" class="mathjax_ignore" alttext="{\displaystyle \phi }"> <semantics> <mrow class="MJX-TeXAtom-ORD"> <mstyle displaystyle="true" scriptlevel="0"> <mi>ϕ<!-- ϕ --></mi> </mstyle> </mrow> <annotation encoding="application/x-tex">{\displaystyle \phi }</annotation> </semantics> </math></span><img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/72b1f30316670aee6270a28334bdf4f5072cdde4" class="mwe-math-fallback-image-inline mw-invert skin-invert" aria-hidden="true" style="vertical-align: -0.671ex; width:1.385ex; height:2.509ex;" alt="{\displaystyle \phi }"/></span> to the solution <b id="mwjQ">repeat</b> until convergence <b id="mwjg">for</b> <span class="texhtml mvar" style="font-style:italic;" about="#mwt55" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"mvar","href":"./Template:Mvar"},"params":{"1":{"wt":"i"}},"i":0}}]}' id="mwjw">i</span> <b id="mwkA">from</b> 1 <b id="mwkQ">until</b> <span class="texhtml mvar" style="font-style:italic;" about="#mwt56" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"mvar","href":"./Template:Mvar"},"params":{"1":{"wt":"n"}},"i":0}}]}' id="mwkg">n</span> <b id="mwkw">do</b> set <span class="texhtml mvar" style="font-style:italic;" about="#mwt57" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"mvar","href":"./Template:Mvar"},"params":{"1":{"wt":"σ"}},"i":0}}]}' id="mwlA">σ</span> to 0 <b id="mwlQ">for</b> <span class="texhtml mvar" style="font-style:italic;" about="#mwt58" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"mvar","href":"./Template:Mvar"},"params":{"1":{"wt":"j"}},"i":0}}]}' id="mwlg">j</span> <b id="mwlw">from</b> 1 <b id="mwmA">until</b> <span class="texhtml mvar" style="font-style:italic;" about="#mwt59" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"mvar","href":"./Template:Mvar"},"params":{"1":{"wt":"n"}},"i":0}}]}' id="mwmQ">n</span> <b id="mwmg">do</b> <b id="mwmw">if</b> <span class="texhtml mvar" style="font-style:italic;" about="#mwt60" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"mvar","href":"./Template:Mvar"},"params":{"1":{"wt":"j"}},"i":0}}]}' id="mwnA">j</span> <span typeof="mw:Entity" id="mwnQ">≠</span> <span class="texhtml mvar" style="font-style:italic;" about="#mwt61" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"mvar","href":"./Template:Mvar"},"params":{"1":{"wt":"i"}},"i":0}}]}' id="mwng">i</span> <b id="mwnw">then</b> set <span class="texhtml mvar" style="font-style:italic;" about="#mwt62" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"mvar","href":"./Template:Mvar"},"params":{"1":{"wt":"σ"}},"i":0}}]}' id="mwoA">σ</span> to <span class="mwe-math-element mwe-math-element-inline" typeof="mw:Extension/math" about="#mwt63" data-mw='{"name":"math","attrs":{},"body":{"extsrc":"\\sigma + a_{ij} \\phi_{j}"}}' id="mwoQ"><span class="mwe-math-mathml-a11y" style="display: none;"><math xmlns="http://www.w3.org/1998/Math/MathML" class="mathjax_ignore" alttext="{\displaystyle \sigma +a_{ij}\phi _{j}}"> <semantics> <mrow class="MJX-TeXAtom-ORD"> <mstyle displaystyle="true" scriptlevel="0"> <mi>σ<!-- σ --></mi> <mo>+</mo> <msub> <mi>a</mi> <mrow class="MJX-TeXAtom-ORD"> <mi>i</mi> <mi>j</mi> </mrow> </msub> <msub> <mi>ϕ<!-- ϕ --></mi> <mrow class="MJX-TeXAtom-ORD"> <mi>j</mi> </mrow> </msub> </mstyle> </mrow> <annotation encoding="application/x-tex">{\displaystyle \sigma +a_{ij}\phi _{j}}</annotation> </semantics> </math></span><img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/ba0617b5c747536cc1046d518c7abd1a6464d963" class="mwe-math-fallback-image-inline mw-invert skin-invert" aria-hidden="true" style="vertical-align: -1.005ex; width:9.172ex; height:2.843ex;" alt="{\displaystyle \sigma +a_{ij}\phi _{j}}"/></span> <b id="mwog">end if</b> <b id="mwow">end</b> (<span class="texhtml mvar" style="font-style:italic;" about="#mwt64" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"mvar","href":"./Template:Mvar"},"params":{"1":{"wt":"j"}},"i":0}}]}' id="mwpA">j</span>-loop) set <span class="mwe-math-element mwe-math-element-inline" typeof="mw:Extension/math" about="#mwt65" data-mw='{"name":"math","attrs":{},"body":{"extsrc":"\\phi_{i}"}}' id="mwpQ"><span class="mwe-math-mathml-a11y" style="display: none;"><math xmlns="http://www.w3.org/1998/Math/MathML" class="mathjax_ignore" alttext="{\displaystyle \phi _{i}}"> <semantics> <mrow class="MJX-TeXAtom-ORD"> <mstyle displaystyle="true" scriptlevel="0"> <msub> <mi>ϕ<!-- ϕ --></mi> <mrow class="MJX-TeXAtom-ORD"> <mi>i</mi> </mrow> </msub> </mstyle> </mrow> <annotation encoding="application/x-tex">{\displaystyle \phi _{i}}</annotation> </semantics> </math></span><img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/0182dbf29b54844c92fd9b0311778a02a38398ec" class="mwe-math-fallback-image-inline mw-invert skin-invert" aria-hidden="true" style="vertical-align: -0.671ex; width:2.185ex; height:2.509ex;" alt="{\displaystyle \phi _{i}}"/></span> to <span class="mwe-math-element mwe-math-element-inline" typeof="mw:Extension/math" about="#mwt66" data-mw='{"name":"math","attrs":{},"body":{"extsrc":"(1-\\omega)\\phi_{i} + \\omega\\,(b_{i}-\\sigma)/a_{ii}"}}' id="mwpg"><span class="mwe-math-mathml-a11y" style="display: none;"><math xmlns="http://www.w3.org/1998/Math/MathML" class="mathjax_ignore" alttext="{\displaystyle (1-\omega )\phi _{i}+\omega \,(b_{i}-\sigma )/a_{ii}}"> <semantics> <mrow class="MJX-TeXAtom-ORD"> <mstyle displaystyle="true" scriptlevel="0"> <mo stretchy="false">(</mo> <mn>1</mn> <mo>−<!-- − --></mo> <mi>ω<!-- ω --></mi> <mo stretchy="false">)</mo> <msub> <mi>ϕ<!-- ϕ --></mi> <mrow class="MJX-TeXAtom-ORD"> <mi>i</mi> </mrow> </msub> <mo>+</mo> <mi>ω<!-- ω --></mi> <mspace width="thinmathspace"></mspace> <mo stretchy="false">(</mo> <msub> <mi>b</mi> <mrow class="MJX-TeXAtom-ORD"> <mi>i</mi> </mrow> </msub> <mo>−<!-- − --></mo> <mi>σ<!-- σ --></mi> <mo stretchy="false">)</mo> <mrow class="MJX-TeXAtom-ORD"> <mo>/</mo> </mrow> <msub> <mi>a</mi> <mrow class="MJX-TeXAtom-ORD"> <mi>i</mi> <mi>i</mi> </mrow> </msub> </mstyle> </mrow> <annotation encoding="application/x-tex">{\displaystyle (1-\omega )\phi _{i}+\omega \,(b_{i}-\sigma )/a_{ii}}</annotation> </semantics> </math></span><img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/d8d4182eacba1122d363934555ed8b855ca4844a" class="mwe-math-fallback-image-inline mw-invert skin-invert" aria-hidden="true" style="vertical-align: -0.838ex; width:25.652ex; height:2.843ex;" alt="{\displaystyle (1-\omega )\phi _{i}+\omega \,(b_{i}-\sigma )/a_{ii}}"/></span> <b id="mwpw">end</b> (<span class="texhtml mvar" style="font-style:italic;" about="#mwt67" typeof="mw:Transclusion" data-mw='{"parts":[{"template":{"target":{"wt":"mvar","href":"./Template:Mvar"},"params":{"1":{"wt":"i"}},"i":0}}]}' id="mwqA">i</span>-loop) check if convergence is reached <b id="mwqQ">end</b> (repeat)

Note: $(1 - ω) ϕ_{i} + \frac{ω}{a_{i i}} (b_{i} - σ)$ can also be written $ϕ_{i} + ω (\frac{b_{i} - σ}{a_{i i}} - ϕ_{i})$ , thus saving one multiplication in each iteration of the outer for-loop.

Example

We are presented the linear system

$\begin{matrix} 4 x_{1} - x_{2} - 6 x_{3} + 0 x_{4} & = 2, \\ - 5 x_{1} - 4 x_{2} + 10 x_{3} + 8 x_{4} & = 21, \\ 0 x_{1} + 9 x_{2} + 4 x_{3} - 2 x_{4} & = - 12, \\ 1 x_{1} + 0 x_{2} - 7 x_{3} + 5 x_{4} & = - 6. \end{matrix}$

To solve the equations, we choose a relaxation factor $ω = 0.5$ and an initial guess vector $ϕ = (0, 0, 0, 0)$ . According to the successive over-relaxation algorithm, the following table is obtained, representing an exemplary iteration with approximations, which ideally, but not necessarily, finds the exact solution, (3, −2, 2, 1), in 38 steps.

A simple implementation of the algorithm in Common Lisp is offered below.

<span id="mw9A"></span><span class="c1" id="mw9Q">;; Set the default floating-point format to "long-float" in order to</span> <span class="c1" id="mw9g">;; ensure correct operation on a wider range of numbers.</span> <span class="p" id="mw9w">(</span><span class="nb" id="mw-A">setf</span><span class="w" id="mw-Q"> </span><span class="vg" id="mw-g">*read-default-float-format*</span><span class="w" id="mw-w"> </span><span class="ss" id="mw_A">'long-float</span><span class="p" id="mw_Q">)</span> <span class="p" id="mw_g">(</span><span class="nb" id="mw_w">defparameter</span><span class="w" id="mwAQA"> </span><span class="nv" id="mwAQE">+MAXIMUM-NUMBER-OF-ITERATIONS+</span><span class="w" id="mwAQI"> </span><span class="mi" id="mwAQM">100</span> <span class="w" id="mwAQQ"> </span><span class="s" id="mwAQU">"The number of iterations beyond which the algorithm should cease its</span> <span class="s" id="mwAQY"> operation, regardless of its current solution. A higher number of</span> <span class="s" id="mwAQc"> iterations might provide a more accurate result, but imposes higher</span> <span class="s" id="mwAQg"> performance requirements."</span><span class="p" id="mwAQk">)</span> <span class="p" id="mwAQo">(</span><span class="nb" id="mwAQs">declaim</span><span class="w" id="mwAQw"> </span><span class="p" id="mwAQ0">(</span><span class="k" id="mwAQ4">type</span><span class="w" id="mwAQ8"> </span><span class="p" id="mwARA">(</span><span class="nc" id="mwARE">integer</span><span class="w" id="mwARI"> </span><span class="mi" id="mwARM">0</span><span class="w" id="mwARQ"> </span><span class="nb" id="mwARU">*</span><span class="p" id="mwARY">)</span><span class="w" id="mwARc"> </span><span class="nv" id="mwARg">+MAXIMUM-NUMBER-OF-ITERATIONS+</span><span class="p" id="mwARk">))</span> <span class="p" id="mwARo">(</span><span class="nb" id="mwARs">defun</span><span class="w" id="mwARw"> </span><span class="nv" id="mwAR0">get-errors</span><span class="w" id="mwAR4"> </span><span class="p" id="mwAR8">(</span><span class="nv" id="mwASA">computed-solution</span><span class="w" id="mwASE"> </span><span class="nv" id="mwASI">exact-solution</span><span class="p" id="mwASM">)</span> <span class="w" id="mwASQ"> </span><span class="s" id="mwASU">"For each component of the COMPUTED-SOLUTION vector, retrieves its</span> <span class="s" id="mwASY"> error with respect to the expected EXACT-SOLUTION vector, returning a</span> <span class="s" id="mwASc"> vector of error values.</span> <span class="s" id="mwASg"> ---</span> <span class="s" id="mwASk"> While both input vectors should be equal in size, this condition is</span> <span class="s" id="mwASo"> not checked and the shortest of the twain determines the output</span> <span class="s" id="mwASs"> vector's number of elements.</span> <span class="s" id="mwASw"> ---</span> <span class="s" id="mwAS0"> The established formula is the following:</span> <span class="s" id="mwAS4"> Let resultVectorSize = min(computedSolution.length, exactSolution.length)</span> <span class="s" id="mwAS8"> Let resultVector = new vector of resultVectorSize</span> <span class="s" id="mwATA"> For i from 0 to (resultVectorSize - 1)</span> <span class="s" id="mwATE"> resultVector[i] = exactSolution[i] - computedSolution[i]</span> <span class="s" id="mwATI"> Return resultVector"</span> <span class="w" id="mwATM"> </span><span class="p" id="mwATQ">(</span><span class="k" id="mwATU">declare</span><span class="w" id="mwATY"> </span><span class="p" id="mwATc">(</span><span class="k" id="mwATg">type</span><span class="w" id="mwATk"> </span><span class="p" id="mwATo">(</span><span class="nb" id="mwATs">vector</span><span class="w" id="mwATw"> </span><span class="nc" id="mwAT0">number</span><span class="w" id="mwAT4"> </span><span class="nb" id="mwAT8">*</span><span class="p" id="mwAUA">)</span><span class="w" id="mwAUE"> </span><span class="nv" id="mwAUI">computed-solution</span><span class="p" id="mwAUM">))</span> <span class="w" id="mwAUQ"> </span><span class="p" id="mwAUU">(</span><span class="k" id="mwAUY">declare</span><span class="w" id="mwAUc"> </span><span class="p" id="mwAUg">(</span><span class="k" id="mwAUk">type</span><span class="w" id="mwAUo"> </span><span class="p" id="mwAUs">(</span><span class="nb" id="mwAUw">vector</span><span class="w" id="mwAU0"> </span><span class="nc" id="mwAU4">number</span><span class="w" id="mwAU8"> </span><span class="nb" id="mwAVA">*</span><span class="p" id="mwAVE">)</span><span class="w" id="mwAVI"> </span><span class="nv" id="mwAVM">exact-solution</span><span class="p" id="mwAVQ">))</span> <span class="w" id="mwAVU"> </span><span class="p" id="mwAVY">(</span><span class="nb" id="mwAVc">map</span><span class="w" id="mwAVg"> </span><span class="o" id="mwAVk">'</span><span class="p" id="mwAVo">(</span><span class="nb" id="mwAVs">vector</span><span class="w" id="mwAVw"> </span><span class="nc" id="mwAV0">number</span><span class="w" id="mwAV4"> </span><span class="nb" id="mwAV8">*</span><span class="p" id="mwAWA">)</span><span class="w" id="mwAWE"> </span><span class="nf" id="mwAWI">#'</span><span class="nb" id="mwAWM">-</span><span class="w" id="mwAWQ"> </span><span class="nv" id="mwAWU">exact-solution</span><span class="w" id="mwAWY"> </span><span class="nv" id="mwAWc">computed-solution</span><span class="p" id="mwAWg">))</span> <span class="p" id="mwAWk">(</span><span class="nb" id="mwAWo">defun</span><span class="w" id="mwAWs"> </span><span class="nv" id="mwAWw">is-convergent</span><span class="w" id="mwAW0"> </span><span class="p" id="mwAW4">(</span><span class="nv" id="mwAW8">errors</span><span class="w" id="mwAXA"> </span><span class="k" id="mwAXE">&key</span><span class="w" id="mwAXI"> </span><span class="p" id="mwAXM">(</span><span class="nv" id="mwAXQ">error-tolerance</span><span class="w" id="mwAXU"> </span><span class="mf" id="mwAXY">0.001</span><span class="p" id="mwAXc">))</span> <span class="w" id="mwAXg"> </span><span class="s" id="mwAXk">"Checks whether the convergence is reached with respect to the</span> <span class="s" id="mwAXo"> ERRORS vector which registers the discrepancy betwixt the computed</span> <span class="s" id="mwAXs"> and the exact solution vector.</span> <span class="s" id="mwAXw"> ---</span> <span class="s" id="mwAX0"> The convergence is fulfilled if and only if each absolute error</span> <span class="s" id="mwAX4"> component is less than or equal to the ERROR-TOLERANCE, that is:</span> <span class="s" id="mwAX8"> For all e in ERRORS, it holds: abs(e) <= errorTolerance."</span> <span class="w" id="mwAYA"> </span><span class="p" id="mwAYE">(</span><span class="k" id="mwAYI">declare</span><span class="w" id="mwAYM"> </span><span class="p" id="mwAYQ">(</span><span class="k" id="mwAYU">type</span><span class="w" id="mwAYY"> </span><span class="p" id="mwAYc">(</span><span class="nb" id="mwAYg">vector</span><span class="w" id="mwAYk"> </span><span class="nc" id="mwAYo">number</span><span class="w" id="mwAYs"> </span><span class="nb" id="mwAYw">*</span><span class="p" id="mwAY0">)</span><span class="w" id="mwAY4"> </span><span class="nv" id="mwAY8">errors</span><span class="p" id="mwAZA">))</span> <span class="w" id="mwAZE"> </span><span class="p" id="mwAZI">(</span><span class="k" id="mwAZM">declare</span><span class="w" id="mwAZQ"> </span><span class="p" id="mwAZU">(</span><span class="k" id="mwAZY">type</span><span class="w" id="mwAZc"> </span><span class="nc" id="mwAZg">number</span><span class="w" id="mwAZk"> </span><span class="nv" id="mwAZo">error-tolerance</span><span class="p" id="mwAZs">))</span> <span class="w" id="mwAZw"> </span><span class="p" id="mwAZ0">(</span><span class="k" id="mwAZ4">flet</span><span class="w" id="mwAZ8"> </span><span class="p" id="mwAaA">((</span><span class="nv" id="mwAaE">error-is-acceptable</span><span class="w" id="mwAaI"> </span><span class="p" id="mwAaM">(</span><span class="nb" id="mwAaQ">error</span><span class="p" id="mwAaU">)</span> <span class="w" id="mwAaY"> </span><span class="p" id="mwAac">(</span><span class="k" id="mwAag">declare</span><span class="w" id="mwAak"> </span><span class="p" id="mwAao">(</span><span class="k" id="mwAas">type</span><span class="w" id="mwAaw"> </span><span class="nc" id="mwAa0">number</span><span class="w" id="mwAa4"> </span><span class="nb" id="mwAa8">error</span><span class="p" id="mwAbA">))</span> <span class="w" id="mwAbE"> </span><span class="p" id="mwAbI">(</span><span class="nb" id="mwAbM"><=</span><span class="w" id="mwAbQ"> </span><span class="p" id="mwAbU">(</span><span class="nb" id="mwAbY">abs</span><span class="w" id="mwAbc"> </span><span class="nb" id="mwAbg">error</span><span class="p" id="mwAbk">)</span><span class="w" id="mwAbo"> </span><span class="nv" id="mwAbs">error-tolerance</span><span class="p" id="mwAbw">)))</span> <span class="w" id="mwAb0"> </span><span class="p" id="mwAb4">(</span><span class="nb" id="mwAb8">every</span><span class="w" id="mwAcA"> </span><span class="nf" id="mwAcE">#'</span><span class="nv" id="mwAcI">error-is-acceptable</span><span class="w" id="mwAcM"> </span><span class="nv" id="mwAcQ">errors</span><span class="p" id="mwAcU">)))</span> <span class="p" id="mwAcY">(</span><span class="nb" id="mwAcc">defun</span><span class="w" id="mwAcg"> </span><span class="nv" id="mwAck">make-zero-vector</span><span class="w" id="mwAco"> </span><span class="p" id="mwAcs">(</span><span class="nv" id="mwAcw">size</span><span class="p" id="mwAc0">)</span> <span class="w" id="mwAc4"> </span><span class="s" id="mwAc8">"Creates and returns a vector of the SIZE with all elements set to 0."</span> <span class="w" id="mwAdA"> </span><span class="p" id="mwAdE">(</span><span class="k" id="mwAdI">declare</span><span class="w" id="mwAdM"> </span><span class="p" id="mwAdQ">(</span><span class="k" id="mwAdU">type</span><span class="w" id="mwAdY"> </span><span class="p" id="mwAdc">(</span><span class="nc" id="mwAdg">integer</span><span class="w" id="mwAdk"> </span><span class="mi" id="mwAdo">0</span><span class="w" id="mwAds"> </span><span class="nb" id="mwAdw">*</span><span class="p" id="mwAd0">)</span><span class="w" id="mwAd4"> </span><span class="nv" id="mwAd8">size</span><span class="p" id="mwAeA">))</span> <span class="w" id="mwAeE"> </span><span class="p" id="mwAeI">(</span><span class="nb" id="mwAeM">make-array</span><span class="w" id="mwAeQ"> </span><span class="nv" id="mwAeU">size</span><span class="w" id="mwAeY"> </span><span class="ss" id="mwAec">:initial-element</span><span class="w" id="mwAeg"> </span><span class="mf" id="mwAek">0.0</span><span class="w" id="mwAeo"> </span><span class="ss" id="mwAes">:element-type</span><span class="w" id="mwAew"> </span><span class="ss" id="mwAe0">'number</span><span class="p" id="mwAe4">))</span> <span class="p" id="mwAe8">(</span><span class="nb" id="mwAfA">defun</span><span class="w" id="mwAfE"> </span><span class="nv" id="mwAfI">successive-over-relaxation</span><span class="w" id="mwAfM"> </span><span class="p" id="mwAfQ">(</span><span class="nv" id="mwAfU">A</span><span class="w" id="mwAfY"> </span><span class="nv" id="mwAfc">b</span><span class="w" id="mwAfg"> </span><span class="nv" id="mwAfk">omega</span> <span class="w" id="mwAfo"> </span><span class="k" id="mwAfs">&key</span><span class="w" id="mwAfw"> </span><span class="p" id="mwAf0">(</span><span class="nv" id="mwAf4">phi</span><span class="w" id="mwAf8"> </span><span class="p" id="mwAgA">(</span><span class="nv" id="mwAgE">make-zero-vector</span><span class="w" id="mwAgI"> </span><span class="p" id="mwAgM">(</span><span class="nb" id="mwAgQ">length</span><span class="w" id="mwAgU"> </span><span class="nv" id="mwAgY">b</span><span class="p" id="mwAgc">)))</span> <span class="w" id="mwAgg"> </span><span class="p" id="mwAgk">(</span><span class="nv" id="mwAgo">convergence-check</span> <span class="w" id="mwAgs"> </span><span class="nf" id="mwAgw">#'</span><span class="p" id="mwAg0">(</span><span class="k" id="mwAg4">lambda</span><span class="w" id="mwAg8"> </span><span class="p" id="mwAhA">(</span><span class="nv" id="mwAhE">iteration</span><span class="w" id="mwAhI"> </span><span class="nv" id="mwAhM">phi</span><span class="p" id="mwAhQ">)</span> <span class="w" id="mwAhU"> </span><span class="p" id="mwAhY">(</span><span class="k" id="mwAhc">declare</span><span class="w" id="mwAhg"> </span><span class="p" id="mwAhk">(</span><span class="k" id="mwAho">ignore</span><span class="w" id="mwAhs"> </span><span class="nv" id="mwAhw">phi</span><span class="p" id="mwAh0">))</span> <span class="w" id="mwAh4"> </span><span class="p" id="mwAh8">(</span><span class="nb" id="mwAiA">>=</span><span class="w" id="mwAiE"> </span><span class="nv" id="mwAiI">iteration</span><span class="w" id="mwAiM"> </span><span class="nv" id="mwAiQ">+MAXIMUM-NUMBER-OF-ITERATIONS+</span><span class="p" id="mwAiU">))))</span> <span class="w" id="mwAiY"> </span><span class="s" id="mwAic">"Implements the successive over-relaxation (SOR) method, applied upon</span> <span class="s" id="mwAig"> the linear equations defined by the matrix A and the right-hand side</span> <span class="s" id="mwAik"> vector B, employing the relaxation factor OMEGA, returning the</span> <span class="s" id="mwAio"> calculated solution vector.</span> <span class="s" id="mwAis"> ---</span> <span class="s" id="mwAiw"> The first algorithm step, the choice of an initial guess PHI, is</span> <span class="s" id="mwAi0"> represented by the optional keyword parameter PHI, which defaults</span> <span class="s" id="mwAi4"> to a zero-vector of the same structure as B. If supplied, this</span> <span class="s" id="mwAi8"> vector will be destructively modified. In any case, the PHI vector</span> <span class="s" id="mwAjA"> constitutes the function's result value.</span> <span class="s" id="mwAjE"> ---</span> <span class="s" id="mwAjI"> The terminating condition is implemented by the CONVERGENCE-CHECK,</span> <span class="s" id="mwAjM"> an optional predicate</span> <span class="s" id="mwAjQ"> lambda(iteration phi) => generalized-boolean</span> <span class="s" id="mwAjU"> which returns T, signifying the immediate termination, upon achieving</span> <span class="s" id="mwAjY"> convergence, or NIL, signaling continuant operation, otherwise. In</span> <span class="s" id="mwAjc"> its default configuration, the CONVERGENCE-CHECK simply abides the</span> <span class="s" id="mwAjg"> iteration's ascension to the ``+MAXIMUM-NUMBER-OF-ITERATIONS+'',</span> <span class="s" id="mwAjk"> ignoring the achieved accuracy of the vector PHI."</span> <span class="w" id="mwAjo"> </span><span class="p" id="mwAjs">(</span><span class="k" id="mwAjw">declare</span><span class="w" id="mwAj0"> </span><span class="p" id="mwAj4">(</span><span class="k" id="mwAj8">type</span><span class="w" id="mwAkA"> </span><span class="p" id="mwAkE">(</span><span class="nc" id="mwAkI">array</span><span class="w" id="mwAkM"> </span><span class="nc" id="mwAkQ">number</span><span class="w" id="mwAkU"> </span><span class="p" id="mwAkY">(</span><span class="nb" id="mwAkc">*</span><span class="w" id="mwAkg"> </span><span class="nb" id="mwAkk">*</span><span class="p" id="mwAko">))</span><span class="w" id="mwAks"> </span><span class="nv" id="mwAkw">A</span><span class="p" id="mwAk0">))</span> <span class="w" id="mwAk4"> </span><span class="p" id="mwAk8">(</span><span class="k" id="mwAlA">declare</span><span class="w" id="mwAlE"> </span><span class="p" id="mwAlI">(</span><span class="k" id="mwAlM">type</span><span class="w" id="mwAlQ"> </span><span class="p" id="mwAlU">(</span><span class="nb" id="mwAlY">vector</span><span class="w" id="mwAlc"> </span><span class="nc" id="mwAlg">number</span><span class="w" id="mwAlk"> </span><span class="nb" id="mwAlo">*</span><span class="p" id="mwAls">)</span><span class="w" id="mwAlw"> </span><span class="nv" id="mwAl0">b</span><span class="p" id="mwAl4">))</span> <span class="w" id="mwAl8"> </span><span class="p" id="mwAmA">(</span><span class="k" id="mwAmE">declare</span><span class="w" id="mwAmI"> </span><span class="p" id="mwAmM">(</span><span class="k" id="mwAmQ">type</span><span class="w" id="mwAmU"> </span><span class="nc" id="mwAmY">number</span><span class="w" id="mwAmc"> </span><span class="nv" id="mwAmg">omega</span><span class="p" id="mwAmk">))</span> <span class="w" id="mwAmo"> </span><span class="p" id="mwAms">(</span><span class="k" id="mwAmw">declare</span><span class="w" id="mwAm0"> </span><span class="p" id="mwAm4">(</span><span class="k" id="mwAm8">type</span><span class="w" id="mwAnA"> </span><span class="p" id="mwAnE">(</span><span class="nb" id="mwAnI">vector</span><span class="w" id="mwAnM"> </span><span class="nc" id="mwAnQ">number</span><span class="w" id="mwAnU"> </span><span class="nb" id="mwAnY">*</span><span class="p" id="mwAnc">)</span><span class="w" id="mwAng"> </span><span class="nv" id="mwAnk">phi</span><span class="p" id="mwAno">))</span> <span class="w" id="mwAns"> </span><span class="p" id="mwAnw">(</span><span class="k" id="mwAn0">declare</span><span class="w" id="mwAn4"> </span><span class="p" id="mwAn8">(</span><span class="k" id="mwAoA">type</span><span class="w" id="mwAoE"> </span><span class="p" id="mwAoI">(</span><span class="k" id="mwAoM">function</span><span class="w" id="mwAoQ"> </span><span class="p" id="mwAoU">((</span><span class="nc" id="mwAoY">integer</span><span class="w" id="mwAoc"> </span><span class="mi" id="mwAog">1</span><span class="w" id="mwAok"> </span><span class="nb" id="mwAoo">*</span><span class="p" id="mwAos">)</span> <span class="w" id="mwAow"> </span><span class="p" id="mwAo0">(</span><span class="nb" id="mwAo4">vector</span><span class="w" id="mwAo8"> </span><span class="nc" id="mwApA">number</span><span class="w" id="mwApE"> </span><span class="nb" id="mwApI">*</span><span class="p" id="mwApM">))</span> <span class="w" id="mwApQ"> </span><span class="nb" id="mwApU">*</span><span class="p" id="mwApY">)</span> <span class="w" id="mwApc"> </span><span class="nv" id="mwApg">convergence-check</span><span class="p" id="mwApk">))</span> <span class="w" id="mwApo"> </span><span class="p" id="mwAps">(</span><span class="k" id="mwApw">let</span><span class="w" id="mwAp0"> </span><span class="p" id="mwAp4">((</span><span class="nv" id="mwAp8">n</span><span class="w" id="mwAqA"> </span><span class="p" id="mwAqE">(</span><span class="nb" id="mwAqI">array-dimension</span><span class="w" id="mwAqM"> </span><span class="nv" id="mwAqQ">A</span><span class="w" id="mwAqU"> </span><span class="mi" id="mwAqY">0</span><span class="p" id="mwAqc">)))</span> <span class="w" id="mwAqg"> </span><span class="p" id="mwAqk">(</span><span class="k" id="mwAqo">declare</span><span class="w" id="mwAqs"> </span><span class="p" id="mwAqw">(</span><span class="k" id="mwAq0">type</span><span class="w" id="mwAq4"> </span><span class="p" id="mwAq8">(</span><span class="nc" id="mwArA">integer</span><span class="w" id="mwArE"> </span><span class="mi" id="mwArI">0</span><span class="w" id="mwArM"> </span><span class="nb" id="mwArQ">*</span><span class="p" id="mwArU">)</span><span class="w" id="mwArY"> </span><span class="nv" id="mwArc">n</span><span class="p" id="mwArg">))</span> <span class="w" id="mwArk"> </span><span class="p" id="mwAro">(</span><span class="nb" id="mwArs">loop</span><span class="w" id="mwArw"> </span><span class="nv" id="mwAr0">for</span><span class="w" id="mwAr4"> </span><span class="nv" id="mwAr8">iteration</span><span class="w" id="mwAsA"> </span><span class="nv" id="mwAsE">from</span><span class="w" id="mwAsI"> </span><span class="mi" id="mwAsM">1</span><span class="w" id="mwAsQ"> </span><span class="nv" id="mwAsU">by</span><span class="w" id="mwAsY"> </span><span class="mi" id="mwAsc">1</span><span class="w" id="mwAsg"> </span><span class="nb" id="mwAsk">do</span> <span class="w" id="mwAso"> </span><span class="p" id="mwAss">(</span><span class="nb" id="mwAsw">loop</span><span class="w" id="mwAs0"> </span><span class="nv" id="mwAs4">for</span><span class="w" id="mwAs8"> </span><span class="nv" id="mwAtA">i</span><span class="w" id="mwAtE"> </span><span class="nv" id="mwAtI">from</span><span class="w" id="mwAtM"> </span><span class="mi" id="mwAtQ">0</span><span class="w" id="mwAtU"> </span><span class="nv" id="mwAtY">below</span><span class="w" id="mwAtc"> </span><span class="nv" id="mwAtg">n</span><span class="w" id="mwAtk"> </span><span class="nv" id="mwAto">by</span><span class="w" id="mwAts"> </span><span class="mi" id="mwAtw">1</span><span class="w" id="mwAt0"> </span><span class="nb" id="mwAt4">do</span> <span class="w" id="mwAt8"> </span><span class="p" id="mwAuA">(</span><span class="k" id="mwAuE">let</span><span class="w" id="mwAuI"> </span><span class="p" id="mwAuM">((</span><span class="nv" id="mwAuQ">rho</span><span class="w" id="mwAuU"> </span><span class="mi" id="mwAuY">0</span><span class="p" id="mwAuc">))</span> <span class="w" id="mwAug"> </span><span class="p" id="mwAuk">(</span><span class="k" id="mwAuo">declare</span><span class="w" id="mwAus"> </span><span class="p" id="mwAuw">(</span><span class="k" id="mwAu0">type</span><span class="w" id="mwAu4"> </span><span class="nc" id="mwAu8">number</span><span class="w" id="mwAvA"> </span><span class="nv" id="mwAvE">rho</span><span class="p" id="mwAvI">))</span> <span class="w" id="mwAvM"> </span><span class="p" id="mwAvQ">(</span><span class="nb" id="mwAvU">loop</span><span class="w" id="mwAvY"> </span><span class="nv" id="mwAvc">for</span><span class="w" id="mwAvg"> </span><span class="nv" id="mwAvk">j</span><span class="w" id="mwAvo"> </span><span class="nv" id="mwAvs">from</span><span class="w" id="mwAvw"> </span><span class="mi" id="mwAv0">0</span><span class="w" id="mwAv4"> </span><span class="nv" id="mwAv8">below</span><span class="w" id="mwAwA"> </span><span class="nv" id="mwAwE">n</span><span class="w" id="mwAwI"> </span><span class="nv" id="mwAwM">by</span><span class="w" id="mwAwQ"> </span><span class="mi" id="mwAwU">1</span><span class="w" id="mwAwY"> </span><span class="nb" id="mwAwc">do</span> <span class="w" id="mwAwg"> </span><span class="p" id="mwAwk">(</span><span class="nb" id="mwAwo">when</span><span class="w" id="mwAws"> </span><span class="p" id="mwAww">(</span><span class="nb" id="mwAw0">/=</span><span class="w" id="mwAw4"> </span><span class="nv" id="mwAw8">j</span><span class="w" id="mwAxA"> </span><span class="nv" id="mwAxE">i</span><span class="p" id="mwAxI">)</span> <span class="w" id="mwAxM"> </span><span class="p" id="mwAxQ">(</span><span class="k" id="mwAxU">let</span><span class="w" id="mwAxY"> </span><span class="p" id="mwAxc">((</span><span class="nv" id="mwAxg">a[ij]</span><span class="w" id="mwAxk"> </span><span class="p" id="mwAxo">(</span><span class="nb" id="mwAxs">aref</span><span class="w" id="mwAxw"> </span><span class="nv" id="mwAx0">A</span><span class="w" id="mwAx4"> </span><span class="nv" id="mwAx8">i</span><span class="w" id="mwAyA"> </span><span class="nv" id="mwAyE">j</span><span class="p" id="mwAyI">))</span> <span class="w" id="mwAyM"> </span><span class="p" id="mwAyQ">(</span><span class="nv" id="mwAyU">phi[j]</span><span class="w" id="mwAyY"> </span><span class="p" id="mwAyc">(</span><span class="nb" id="mwAyg">aref</span><span class="w" id="mwAyk"> </span><span class="nv" id="mwAyo">phi</span><span class="w" id="mwAys"> </span><span class="nv" id="mwAyw">j</span><span class="p" id="mwAy0">)))</span> <span class="w" id="mwAy4"> </span><span class="p" id="mwAy8">(</span><span class="nb" id="mwAzA">incf</span><span class="w" id="mwAzE"> </span><span class="nv" id="mwAzI">rho</span><span class="w" id="mwAzM"> </span><span class="p" id="mwAzQ">(</span><span class="nb" id="mwAzU">*</span><span class="w" id="mwAzY"> </span><span class="nv" id="mwAzc">a[ij]</span><span class="w" id="mwAzg"> </span><span class="nv" id="mwAzk">phi[j]</span><span class="p" id="mwAzo">)))))</span> <span class="w" id="mwAzs"> </span><span class="p" id="mwAzw">(</span><span class="nb" id="mwAz0">setf</span><span class="w" id="mwAz4"> </span><span class="p" id="mwAz8">(</span><span class="nb" id="mwA0A">aref</span><span class="w" id="mwA0E"> </span><span class="nv" id="mwA0I">phi</span><span class="w" id="mwA0M"> </span><span class="nv" id="mwA0Q">i</span><span class="p" id="mwA0U">)</span> <span class="w" id="mwA0Y"> </span><span class="p" id="mwA0c">(</span><span class="nb" id="mwA0g">+</span><span class="w" id="mwA0k"> </span><span class="p" id="mwA0o">(</span><span class="nb" id="mwA0s">*</span><span class="w" id="mwA0w"> </span><span class="p" id="mwA00">(</span><span class="nb" id="mwA04">-</span><span class="w" id="mwA08"> </span><span class="mi" id="mwA1A">1</span><span class="w" id="mwA1E"> </span><span class="nv" id="mwA1I">omega</span><span class="p" id="mwA1M">)</span> <span class="w" id="mwA1Q"> </span><span class="p" id="mwA1U">(</span><span class="nb" id="mwA1Y">aref</span><span class="w" id="mwA1c"> </span><span class="nv" id="mwA1g">phi</span><span class="w" id="mwA1k"> </span><span class="nv" id="mwA1o">i</span><span class="p" id="mwA1s">))</span> <span class="w" id="mwA1w"> </span><span class="p" id="mwA10">(</span><span class="nb" id="mwA14">*</span><span class="w" id="mwA18"> </span><span class="p" id="mwA2A">(</span><span class="nb" id="mwA2E">/</span><span class="w" id="mwA2I"> </span><span class="nv" id="mwA2M">omega</span><span class="w" id="mwA2Q"> </span><span class="p" id="mwA2U">(</span><span class="nb" id="mwA2Y">aref</span><span class="w" id="mwA2c"> </span><span class="nv" id="mwA2g">A</span><span class="w" id="mwA2k"> </span><span class="nv" id="mwA2o">i</span><span class="w" id="mwA2s"> </span><span class="nv" id="mwA2w">i</span><span class="p" id="mwA20">))</span> <span class="w" id="mwA24"> </span><span class="p" id="mwA28">(</span><span class="nb" id="mwA3A">-</span><span class="w" id="mwA3E"> </span><span class="p" id="mwA3I">(</span><span class="nb" id="mwA3M">aref</span><span class="w" id="mwA3Q"> </span><span class="nv" id="mwA3U">b</span><span class="w" id="mwA3Y"> </span><span class="nv" id="mwA3c">i</span><span class="p" id="mwA3g">)</span><span class="w" id="mwA3k"> </span><span class="nv" id="mwA3o">rho</span><span class="p" id="mwA3s">))))))</span> <span class="w" id="mwA3w"> </span><span class="p" id="mwA30">(</span><span class="nb" id="mwA34">format</span><span class="w" id="mwA38"> </span><span class="no" id="mwA4A">T</span><span class="w" id="mwA4E"> </span><span class="s" id="mwA4I">"~&~d. solution = ~a"</span><span class="w" id="mwA4M"> </span><span class="nv" id="mwA4Q">iteration</span><span class="w" id="mwA4U"> </span><span class="nv" id="mwA4Y">phi</span><span class="p" id="mwA4c">)</span> <span class="w" id="mwA4g"> </span><span class="c1" id="mwA4k">;; Check if convergence is reached.</span> <span class="w" id="mwA4o"> </span><span class="p" id="mwA4s">(</span><span class="nb" id="mwA4w">when</span><span class="w" id="mwA40"> </span><span class="p" id="mwA44">(</span><span class="nb" id="mwA48">funcall</span><span class="w" id="mwA5A"> </span><span class="nv" id="mwA5E">convergence-check</span><span class="w" id="mwA5I"> </span><span class="nv" id="mwA5M">iteration</span><span class="w" id="mwA5Q"> </span><span class="nv" id="mwA5U">phi</span><span class="p" id="mwA5Y">)</span> <span class="w" id="mwA5c"> </span><span class="p" id="mwA5g">(</span><span class="nb" id="mwA5k">return</span><span class="p" id="mwA5o">))))</span> <span class="w" id="mwA5s"> </span><span class="p" id="mwA5w">(</span><span class="k" id="mwA50">the</span><span class="w" id="mwA54"> </span><span class="p" id="mwA58">(</span><span class="nb" id="mwA6A">vector</span><span class="w" id="mwA6E"> </span><span class="nc" id="mwA6I">number</span><span class="w" id="mwA6M"> </span><span class="nb" id="mwA6Q">*</span><span class="p" id="mwA6U">)</span><span class="w" id="mwA6Y"> </span><span class="nv" id="mwA6c">phi</span><span class="p" id="mwA6g">))</span> <span class="c1" id="mwA6k">;; Summon the function with the exemplary parameters.</span> <span class="p" id="mwA6o">(</span><span class="k" id="mwA6s">let</span><span class="w" id="mwA6w"> </span><span class="p" id="mwA60">((</span><span class="nv" id="mwA64">A</span><span class="w" id="mwA68"> </span><span class="p" id="mwA7A">(</span><span class="nb" id="mwA7E">make-array</span><span class="w" id="mwA7I"> </span><span class="p" id="mwA7M">(</span><span class="nb" id="mwA7Q">list</span><span class="w" id="mwA7U"> </span><span class="mi" id="mwA7Y">4</span><span class="w" id="mwA7c"> </span><span class="mi" id="mwA7g">4</span><span class="p" id="mwA7k">)</span> <span class="w" id="mwA7o"> </span><span class="ss" id="mwA7s">:initial-contents</span> <span class="w" id="mwA7w"> </span><span class="o" id="mwA70">'</span><span class="p" id="mwA74">((</span><span class="w" id="mwA78"> </span><span class="mi" id="mwA8A">4</span><span class="w" id="mwA8E"> </span><span class="mi" id="mwA8I">-1</span><span class="w" id="mwA8M"> </span><span class="mi" id="mwA8Q">-6</span><span class="w" id="mwA8U"> </span><span class="mi" id="mwA8Y">0</span><span class="w" id="mwA8c"> </span><span class="p" id="mwA8g">)</span> <span class="w" id="mwA8k"> </span><span class="p" id="mwA8o">(</span><span class="w" id="mwA8s"> </span><span class="mi" id="mwA8w">-5</span><span class="w" id="mwA80"> </span><span class="mi" id="mwA84">-4</span><span class="w" id="mwA88"> </span><span class="mi" id="mwA9A">10</span><span class="w" id="mwA9E"> </span><span class="mi" id="mwA9I">8</span><span class="w" id="mwA9M"> </span><span class="p" id="mwA9Q">)</span> <span class="w" id="mwA9U"> </span><span class="p" id="mwA9Y">(</span><span class="w" id="mwA9c"> </span><span class="mi" id="mwA9g">0</span><span class="w" id="mwA9k"> </span><span class="mi" id="mwA9o">9</span><span class="w" id="mwA9s"> </span><span class="mi" id="mwA9w">4</span><span class="w" id="mwA90"> </span><span class="mi" id="mwA94">-2</span><span class="w" id="mwA98"> </span><span class="p" id="mwA-A">)</span> <span class="w" id="mwA-E"> </span><span class="p" id="mwA-I">(</span><span class="w" id="mwA-M"> </span><span class="mi" id="mwA-Q">1</span><span class="w" id="mwA-U"> </span><span class="mi" id="mwA-Y">0</span><span class="w" id="mwA-c"> </span><span class="mi" id="mwA-g">-7</span><span class="w" id="mwA-k"> </span><span class="mi" id="mwA-o">5</span><span class="w" id="mwA-s"> </span><span class="p" id="mwA-w">))))</span> <span class="w" id="mwA-0"> </span><span class="p" id="mwA-4">(</span><span class="nv" id="mwA-8">b</span><span class="w" id="mwA_A"> </span><span class="p" id="mwA_E">(</span><span class="nb" id="mwA_I">vector</span><span class="w" id="mwA_M"> </span><span class="mi" id="mwA_Q">2</span><span class="w" id="mwA_U"> </span><span class="mi" id="mwA_Y">21</span><span class="w" id="mwA_c"> </span><span class="mi" id="mwA_g">-12</span><span class="w" id="mwA_k"> </span><span class="mi" id="mwA_o">-6</span><span class="p" id="mwA_s">))</span> <span class="w" id="mwA_w"> </span><span class="p" id="mwA_0">(</span><span class="nv" id="mwA_4">omega</span><span class="w" id="mwA_8"> </span><span class="mf" id="mwBAA">0.5</span><span class="p" id="mwBAE">)</span> <span class="w" id="mwBAI"> </span><span class="p" id="mwBAM">(</span><span class="nv" id="mwBAQ">exact-solution</span><span class="w" id="mwBAU"> </span><span class="p" id="mwBAY">(</span><span class="nb" id="mwBAc">vector</span><span class="w" id="mwBAg"> </span><span class="mi" id="mwBAk">3</span><span class="w" id="mwBAo"> </span><span class="mi" id="mwBAs">-2</span><span class="w" id="mwBAw"> </span><span class="mi" id="mwBA0">2</span><span class="w" id="mwBA4"> </span><span class="mi" id="mwBA8">1</span><span class="p" id="mwBBA">)))</span> <span class="w" id="mwBBE"> </span><span class="p" id="mwBBI">(</span><span class="nv" id="mwBBM">successive-over-relaxation</span> <span class="w" id="mwBBQ"> </span><span class="nv" id="mwBBU">A</span><span class="w" id="mwBBY"> </span><span class="nv" id="mwBBc">b</span><span class="w" id="mwBBg"> </span><span class="nv" id="mwBBk">omega</span> <span class="w" id="mwBBo"> </span><span class="ss" id="mwBBs">:convergence-check</span> <span class="w" id="mwBBw"> </span><span class="nf" id="mwBB0">#'</span><span class="p" id="mwBB4">(</span><span class="k" id="mwBB8">lambda</span><span class="w" id="mwBCA"> </span><span class="p" id="mwBCE">(</span><span class="nv" id="mwBCI">iteration</span><span class="w" id="mwBCM"> </span><span class="nv" id="mwBCQ">phi</span><span class="p" id="mwBCU">)</span> <span class="w" id="mwBCY"> </span><span class="p" id="mwBCc">(</span><span class="k" id="mwBCg">declare</span><span class="w" id="mwBCk"> </span><span class="p" id="mwBCo">(</span><span class="k" id="mwBCs">type</span><span class="w" id="mwBCw"> </span><span class="p" id="mwBC0">(</span><span class="nc" id="mwBC4">integer</span><span class="w" id="mwBC8"> </span><span class="mi" id="mwBDA">0</span><span class="w" id="mwBDE"> </span><span class="nb" id="mwBDI">*</span><span class="p" id="mwBDM">)</span><span class="w" id="mwBDQ"> </span><span class="nv" id="mwBDU">iteration</span><span class="p" id="mwBDY">))</span> <span class="w" id="mwBDc"> </span><span class="p" id="mwBDg">(</span><span class="k" id="mwBDk">declare</span><span class="w" id="mwBDo"> </span><span class="p" id="mwBDs">(</span><span class="k" id="mwBDw">type</span><span class="w" id="mwBD0"> </span><span class="p" id="mwBD4">(</span><span class="nb" id="mwBD8">vector</span><span class="w" id="mwBEA"> </span><span class="nc" id="mwBEE">number</span><span class="w" id="mwBEI"> </span><span class="nb" id="mwBEM">*</span><span class="p" id="mwBEQ">)</span><span class="w" id="mwBEU"> </span><span class="nv" id="mwBEY">phi</span><span class="p" id="mwBEc">))</span> <span class="w" id="mwBEg"> </span><span class="p" id="mwBEk">(</span><span class="k" id="mwBEo">let</span><span class="w" id="mwBEs"> </span><span class="p" id="mwBEw">((</span><span class="nv" id="mwBE0">errors</span><span class="w" id="mwBE4"> </span><span class="p" id="mwBE8">(</span><span class="nv" id="mwBFA">get-errors</span><span class="w" id="mwBFE"> </span><span class="nv" id="mwBFI">phi</span><span class="w" id="mwBFM"> </span><span class="nv" id="mwBFQ">exact-solution</span><span class="p" id="mwBFU">)))</span> <span class="w" id="mwBFY"> </span><span class="p" id="mwBFc">(</span><span class="k" id="mwBFg">declare</span><span class="w" id="mwBFk"> </span><span class="p" id="mwBFo">(</span><span class="k" id="mwBFs">type</span><span class="w" id="mwBFw"> </span><span class="p" id="mwBF0">(</span><span class="nb" id="mwBF4">vector</span><span class="w" id="mwBF8"> </span><span class="nc" id="mwBGA">number</span><span class="w" id="mwBGE"> </span><span class="nb" id="mwBGI">*</span><span class="p" id="mwBGM">)</span><span class="w" id="mwBGQ"> </span><span class="nv" id="mwBGU">errors</span><span class="p" id="mwBGY">))</span> <span class="w" id="mwBGc"> </span><span class="p" id="mwBGg">(</span><span class="nb" id="mwBGk">format</span><span class="w" id="mwBGo"> </span><span class="no" id="mwBGs">T</span><span class="w" id="mwBGw"> </span><span class="s" id="mwBG0">"~&~d. errors = ~a"</span><span class="w" id="mwBG4"> </span><span class="nv" id="mwBG8">iteration</span><span class="w" id="mwBHA"> </span><span class="nv" id="mwBHE">errors</span><span class="p" id="mwBHI">)</span> <span class="w" id="mwBHM"> </span><span class="p" id="mwBHQ">(</span><span class="nb" id="mwBHU">or</span><span class="w" id="mwBHY"> </span><span class="p" id="mwBHc">(</span><span class="nv" id="mwBHg">is-convergent</span><span class="w" id="mwBHk"> </span><span class="nv" id="mwBHo">errors</span><span class="w" id="mwBHs"> </span><span class="ss" id="mwBHw">:error-tolerance</span><span class="w" id="mwBH0"> </span><span class="mf" id="mwBH4">0.0</span><span class="p" id="mwBH8">)</span> <span class="w" id="mwBIA"> </span><span class="p" id="mwBIE">(</span><span class="nb" id="mwBII">>=</span><span class="w" id="mwBIM"> </span><span class="nv" id="mwBIQ">iteration</span><span class="w" id="mwBIU"> </span><span class="nv" id="mwBIY">+MAXIMUM-NUMBER-OF-ITERATIONS+</span><span class="p" id="mwBIc">))))))</span>

A simple Python implementation of the pseudo-code provided above.

<span id="mwBIs"></span><span class="kn" id="mwBIw">import</span><span class="w" id="mwBI0"> </span><span class="nn" id="mwBI4">numpy</span><span class="w" id="mwBI8"> </span><span class="k" id="mwBJA">as</span><span class="w" id="mwBJE"> </span><span class="nn" id="mwBJI">np</span> <span class="kn" id="mwBJM">from</span><span class="w" id="mwBJQ"> </span><span class="nn" id="mwBJU">scipy</span><span class="w" id="mwBJY"> </span><span class="kn" id="mwBJc">import</span> <span class="n" id="mwBJg">linalg</span> <span class="k" id="mwBJk">def</span><span class="w" id="mwBJo"> </span><span class="nf" id="mwBJs">sor_solver</span><span class="p" id="mwBJw">(</span><span class="n" id="mwBJ0">A</span><span class="p" id="mwBJ4">,</span> <span class="n" id="mwBJ8">b</span><span class="p" id="mwBKA">,</span> <span class="n" id="mwBKE">omega</span><span class="p" id="mwBKI">,</span> <span class="n" id="mwBKM">initial_guess</span><span class="p" id="mwBKQ">,</span> <span class="n" id="mwBKU">convergence_criteria</span><span class="p" id="mwBKY">):</span> <span class="w" id="mwBKc"> </span><span class="sd" id="mwBKg">"""</span> <span class="sd" id="mwBKk"> This is an implementation of the pseudo-code provided in the Wikipedia article.</span> <span class="sd" id="mwBKo"> Arguments:</span> <span class="sd" id="mwBKs"> A: nxn numpy matrix.</span> <span class="sd" id="mwBKw"> b: n dimensional numpy vector.</span> <span class="sd" id="mwBK0"> omega: relaxation factor.</span> <span class="sd" id="mwBK4"> initial_guess: An initial solution guess for the solver to start with.</span> <span class="sd" id="mwBK8"> convergence_criteria: The maximum discrepancy acceptable to regard the current solution as fitting.</span> <span class="sd" id="mwBLA"> Returns:</span> <span class="sd" id="mwBLE"> phi: solution vector of dimension n.</span> <span class="sd" id="mwBLI"> """</span> <span class="n" id="mwBLM">step</span> <span class="o" id="mwBLQ">=</span> <span class="mi" id="mwBLU">0</span> <span class="n" id="mwBLY">phi</span> <span class="o" id="mwBLc">=</span> <span class="n" id="mwBLg">initial_guess</span><span class="p" id="mwBLk">[:]</span> <span class="n" id="mwBLo">residual</span> <span class="o" id="mwBLs">=</span> <span class="n" id="mwBLw">linalg</span><span class="o" id="mwBL0">.</span><span class="n" id="mwBL4">norm</span><span class="p" id="mwBL8">(</span><span class="n" id="mwBMA">A</span> <span class="o" id="mwBME">@</span> <span class="n" id="mwBMI">phi</span> <span class="o" id="mwBMM">-</span> <span class="n" id="mwBMQ">b</span><span class="p" id="mwBMU">)</span> <span class="c1" id="mwBMY"># Initial residual</span> <span class="k" id="mwBMc">while</span> <span class="n" id="mwBMg">residual</span> <span class="o" id="mwBMk">></span> <span class="n" id="mwBMo">convergence_criteria</span><span class="p" id="mwBMs">:</span> <span class="k" id="mwBMw">for</span> <span class="n" id="mwBM0">i</span> <span class="ow" id="mwBM4">in</span> <span class="nb" id="mwBM8">range</span><span class="p" id="mwBNA">(</span><span class="n" id="mwBNE">A</span><span class="o" id="mwBNI">.</span><span class="n" id="mwBNM">shape</span><span class="p" id="mwBNQ">[</span><span class="mi" id="mwBNU">0</span><span class="p" id="mwBNY">]):</span> <span class="n" id="mwBNc">sigma</span> <span class="o" id="mwBNg">=</span> <span class="mi" id="mwBNk">0</span> <span class="k" id="mwBNo">for</span> <span class="n" id="mwBNs">j</span> <span class="ow" id="mwBNw">in</span> <span class="nb" id="mwBN0">range</span><span class="p" id="mwBN4">(</span><span class="n" id="mwBN8">A</span><span class="o" id="mwBOA">.</span><span class="n" id="mwBOE">shape</span><span class="p" id="mwBOI">[</span><span class="mi" id="mwBOM">1</span><span class="p" id="mwBOQ">]):</span> <span class="k" id="mwBOU">if</span> <span class="n" id="mwBOY">j</span> <span class="o" id="mwBOc">!=</span> <span class="n" id="mwBOg">i</span><span class="p" id="mwBOk">:</span> <span class="n" id="mwBOo">sigma</span> <span class="o" id="mwBOs">+=</span> <span class="n" id="mwBOw">A</span><span class="p" id="mwBO0">[</span><span class="n" id="mwBO4">i</span><span class="p" id="mwBO8">,</span> <span class="n" id="mwBPA">j</span><span class="p" id="mwBPE">]</span> <span class="o" id="mwBPI">*</span> <span class="n" id="mwBPM">phi</span><span class="p" id="mwBPQ">[</span><span class="n" id="mwBPU">j</span><span class="p" id="mwBPY">]</span> <span class="n" id="mwBPc">phi</span><span class="p" id="mwBPg">[</span><span class="n" id="mwBPk">i</span><span class="p" id="mwBPo">]</span> <span class="o" id="mwBPs">=</span> <span class="p" id="mwBPw">(</span><span class="mi" id="mwBP0">1</span> <span class="o" id="mwBP4">-</span> <span class="n" id="mwBP8">omega</span><span class="p" id="mwBQA">)</span> <span class="o" id="mwBQE">*</span> <span class="n" id="mwBQI">phi</span><span class="p" id="mwBQM">[</span><span class="n" id="mwBQQ">i</span><span class="p" id="mwBQU">]</span> <span class="o" id="mwBQY">+</span> <span class="p" id="mwBQc">(</span><span class="n" id="mwBQg">omega</span> <span class="o" id="mwBQk">/</span> <span class="n" id="mwBQo">A</span><span class="p" id="mwBQs">[</span><span class="n" id="mwBQw">i</span><span class="p" id="mwBQ0">,</span> <span class="n" id="mwBQ4">i</span><span class="p" id="mwBQ8">])</span> <span class="o" id="mwBRA">*</span> <span class="p" id="mwBRE">(</span><span class="n" id="mwBRI">b</span><span class="p" id="mwBRM">[</span><span class="n" id="mwBRQ">i</span><span class="p" id="mwBRU">]</span> <span class="o" id="mwBRY">-</span> <span class="n" id="mwBRc">sigma</span><span class="p" id="mwBRg">)</span> <span class="n" id="mwBRk">residual</span> <span class="o" id="mwBRo">=</span> <span class="n" id="mwBRs">linalg</span><span class="o" id="mwBRw">.</span><span class="n" id="mwBR0">norm</span><span class="p" id="mwBR4">(</span><span class="n" id="mwBR8">A</span> <span class="o" id="mwBSA">@</span> <span class="n" id="mwBSE">phi</span> <span class="o" id="mwBSI">-</span> <span class="n" id="mwBSM">b</span><span class="p" id="mwBSQ">)</span> <span class="n" id="mwBSU">step</span> <span class="o" id="mwBSY">+=</span> <span class="mi" id="mwBSc">1</span> <span class="nb" id="mwBSg">print</span><span class="p" id="mwBSk">(</span><span class="s2" id="mwBSo">"Step </span><span class="si" id="mwBSs">{}</span><span class="s2" id="mwBSw"> Residual: </span><span class="si" id="mwBS0">{:10.6g}</span><span class="s2" id="mwBS4">"</span><span class="o" id="mwBS8">.</span><span class="n" id="mwBTA">format</span><span class="p" id="mwBTE">(</span><span class="n" id="mwBTI">step</span><span class="p" id="mwBTM">,</span> <span class="n" id="mwBTQ">residual</span><span class="p" id="mwBTU">))</span> <span class="k" id="mwBTY">return</span> <span class="n" id="mwBTc">phi</span> <span class="c1" id="mwBTg"># An example case that mirrors the one in the Wikipedia article</span> <span class="n" id="mwBTk">residual_convergence</span> <span class="o" id="mwBTo">=</span> <span class="mf" id="mwBTs">1e-8</span> <span class="n" id="mwBTw">omega</span> <span class="o" id="mwBT0">=</span> <span class="mf" id="mwBT4">0.5</span> <span class="c1" id="mwBT8"># Relaxation factor</span> <span class="n" id="mwBUA">A</span> <span class="o" id="mwBUE">=</span> <span class="n" id="mwBUI">np</span><span class="o" id="mwBUM">.</span><span class="n" id="mwBUQ">array</span><span class="p" id="mwBUU">([[</span><span class="mi" id="mwBUY">4</span><span class="p" id="mwBUc">,</span> <span class="o" id="mwBUg">-</span><span class="mi" id="mwBUk">1</span><span class="p" id="mwBUo">,</span> <span class="o" id="mwBUs">-</span><span class="mi" id="mwBUw">6</span><span class="p" id="mwBU0">,</span> <span class="mi" id="mwBU4">0</span><span class="p" id="mwBU8">],</span> <span class="p" id="mwBVA">[</span><span class="o" id="mwBVE">-</span><span class="mi" id="mwBVI">5</span><span class="p" id="mwBVM">,</span> <span class="o" id="mwBVQ">-</span><span class="mi" id="mwBVU">4</span><span class="p" id="mwBVY">,</span> <span class="mi" id="mwBVc">10</span><span class="p" id="mwBVg">,</span> <span class="mi" id="mwBVk">8</span><span class="p" id="mwBVo">],</span> <span class="p" id="mwBVs">[</span><span class="mi" id="mwBVw">0</span><span class="p" id="mwBV0">,</span> <span class="mi" id="mwBV4">9</span><span class="p" id="mwBV8">,</span> <span class="mi" id="mwBWA">4</span><span class="p" id="mwBWE">,</span> <span class="o" id="mwBWI">-</span><span class="mi" id="mwBWM">2</span><span class="p" id="mwBWQ">],</span> <span class="p" id="mwBWU">[</span><span class="mi" id="mwBWY">1</span><span class="p" id="mwBWc">,</span> <span class="mi" id="mwBWg">0</span><span class="p" id="mwBWk">,</span> <span class="o" id="mwBWo">-</span><span class="mi" id="mwBWs">7</span><span class="p" id="mwBWw">,</span> <span class="mi" id="mwBW0">5</span><span class="p" id="mwBW4">]])</span> <span class="n" id="mwBW8">b</span> <span class="o" id="mwBXA">=</span> <span class="n" id="mwBXE">np</span><span class="o" id="mwBXI">.</span><span class="n" id="mwBXM">array</span><span class="p" id="mwBXQ">([</span><span class="mi" id="mwBXU">2</span><span class="p" id="mwBXY">,</span> <span class="mi" id="mwBXc">21</span><span class="p" id="mwBXg">,</span> <span class="o" id="mwBXk">-</span><span class="mi" id="mwBXo">12</span><span class="p" id="mwBXs">,</span> <span class="o" id="mwBXw">-</span><span class="mi" id="mwBX0">6</span><span class="p" id="mwBX4">])</span> <span class="n" id="mwBX8">initial_guess</span> <span class="o" id="mwBYA">=</span> <span class="n" id="mwBYE">np</span><span class="o" id="mwBYI">.</span><span class="n" id="mwBYM">zeros</span><span class="p" id="mwBYQ">(</span><span class="mi" id="mwBYU">4</span><span class="p" id="mwBYY">)</span> <span class="n" id="mwBYc">phi</span> <span class="o" id="mwBYg">=</span> <span class="n" id="mwBYk">sor_solver</span><span class="p" id="mwBYo">(</span><span class="n" id="mwBYs">A</span><span class="p" id="mwBYw">,</span> <span class="n" id="mwBY0">b</span><span class="p" id="mwBY4">,</span> <span class="n" id="mwBY8">omega</span><span class="p" id="mwBZA">,</span> <span class="n" id="mwBZE">initial_guess</span><span class="p" id="mwBZI">,</span> <span class="n" id="mwBZM">residual_convergence</span><span class="p" id="mwBZQ">)</span> <span class="nb" id="mwBZU">print</span><span class="p" id="mwBZY">(</span><span class="n" id="mwBZc">phi</span><span class="p" id="mwBZg">)</span>

Symmetric successive over-relaxation

The version for symmetric matrices A, in which

$U = L^{T},$

is referred to as Symmetric Successive Over-Relaxation, or (SSOR), in which

$P = (\frac{D}{ω} + L) \frac{ω}{2 - ω} D^{- 1} (\frac{D}{ω} + U),$

and the iterative method is

$x^{k + 1} = x^{k} - γ^{k} P^{- 1} (A x^{k} - b), k \geq 0.$

The SOR and SSOR methods are credited to David M. Young Jr.

Other applications of the method

A similar technique can be used for any iterative method. If the original iteration had the form

$x_{n + 1} = f (x_{n})$

then the modified version would use

$x_{n + 1}^{S O R} = (1 - ω) x_{n}^{S O R} + ω f (x_{n}^{S O R}) .$

However, the formulation presented above, used for solving systems of linear equations, is not a special case of this formulation if x is considered to be the complete vector. If this formulation is used instead, the equation for calculating the next vector will look like

$x^{(k + 1)} = (1 - ω) x^{(k)} + ω L_{*}^{- 1} (b - U x^{(k)}),$

where $L_{*} = L + D$ . Values of $ω > 1$ are used to speed up convergence of a slow-converging process, while values of $ω < 1$ are often used to help establish convergence of a diverging iterative process or speed up the convergence of an overshooting process.

There are various methods that adaptively set the relaxation parameter $ω$ based on the observed behavior of the converging process. Usually they help to reach a super-linear convergence for some problems but fail for the others.