.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/linalg/iterative/plot_02_gauss_seidel_and_sor.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_api_gallery_linalg_iterative_plot_02_gauss_seidel_and_sor.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_api_gallery_linalg_iterative_plot_02_gauss_seidel_and_sor.py:


Jacobi, Gauss-Seidel, and successive over-relaxation
=========================================================

The three classical stationary iterations on the 1-D Poisson matrix
``tridiag(-1, 2, -1)``. Gauss-Seidel contracts the error about twice as
fast per sweep as Jacobi, and SOR at Young's optimal relaxation factor is
faster by an order of magnitude.

.. GENERATED FROM PYTHON SOURCE LINES 12-18

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from mathematicskit.linalg import SOR, GaussSeidel, JacobiIteration, jacobi_spectral_radius, optimal_sor_omega
    from mathematicskit.linalg.visualizers.plots import plot_residual_history








.. GENERATED FROM PYTHON SOURCE LINES 19-21

Spectral radii predicted by Young's theory
-----------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 21-31

.. code-block:: Python


    n = 40
    A = 2.0 * np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1)
    b = np.ones(n)
    rho_j = jacobi_spectral_radius(A)
    omega = optimal_sor_omega(A)
    print(f"rho(Jacobi) = {rho_j:.6f}   (cos(pi/(n+1)) = {np.cos(np.pi / (n + 1)):.6f})")
    print(f"rho(Gauss-Seidel) = rho(Jacobi)^2 = {rho_j**2:.6f}")
    print(f"optimal omega = {omega:.6f}, rho(SOR) = omega - 1 = {omega - 1:.6f}")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    rho(Jacobi) = 0.997066   (cos(pi/(n+1)) = 0.997066)
    rho(Gauss-Seidel) = rho(Jacobi)^2 = 0.994140
    optimal omega = 1.857788, rho(SOR) = omega - 1 = 0.857788




.. GENERATED FROM PYTHON SOURCE LINES 32-34

Convergence histories
---------------------------

.. GENERATED FROM PYTHON SOURCE LINES 34-50

.. code-block:: Python


    results = [
        JacobiIteration(tol=1e-8, max_iter=20000).solve(A, b),
        GaussSeidel(tol=1e-8, max_iter=20000).solve(A, b),
        SOR(omega=omega, tol=1e-8, max_iter=20000).solve(A, b),
    ]
    for r in results:
        print(f"{r.method:13s}: {r.iterations:5d} sweeps")

    fig, ax = plt.subplots(figsize=(6, 4))
    for r in results:
        plot_residual_history(r, ax=ax)
    ax.legend()
    fig.tight_layout()

    plt.show()



.. image-sg:: /api/gallery/linalg/iterative/images/sphx_glr_plot_02_gauss_seidel_and_sor_001.png
   :alt: Iterative solver convergence
   :srcset: /api/gallery/linalg/iterative/images/sphx_glr_plot_02_gauss_seidel_and_sor_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    jacobi       :  6238 sweeps
    gauss_seidel :  3120 sweeps
    sor          :   152 sweeps





.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.303 seconds)


.. _sphx_glr_download_api_gallery_linalg_iterative_plot_02_gauss_seidel_and_sor.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_02_gauss_seidel_and_sor.ipynb <plot_02_gauss_seidel_and_sor.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_02_gauss_seidel_and_sor.py <plot_02_gauss_seidel_and_sor.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_02_gauss_seidel_and_sor.zip <plot_02_gauss_seidel_and_sor.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_