.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/pde/elliptic/plot_04_multigrid.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_pde_elliptic_plot_04_multigrid.py>`
        to download the full example code or to run this example in your browser via JupyterLite.

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

.. _sphx_glr_api_gallery_pde_elliptic_plot_04_multigrid.py:


Multigrid: a solver whose cost does not grow with the grid
==========================================================

Relaxation kills wiggly error fast but smooth error slowly, so Jacobi's
iteration count grows like the number of grid points per side squared.
Fedorenko (1964) and Brandt (1977) noticed that smooth error on a fine
grid looks wiggly on a coarser one, where it is cheap to remove. A
multigrid V-cycle smooths, corrects on successively coarser grids, and
smooths again. The number of cycles needed stays flat as the grid is
refined. This script shows Jacobi smoothing a random error, then
compares cycle counts across grid sizes.

.. GENERATED FROM PYTHON SOURCE LINES 16-22

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from mathematicskit.linalg.visualizers import plot_residual_history
    from mathematicskit.pde import multigrid_poisson_2d, solve_poisson_2d








.. GENERATED FROM PYTHON SOURCE LINES 23-27

Relaxation smooths the error
----------------------------
Solving Laplace's equation (exact solution zero) from a random start,
a few damped-Jacobi sweeps leave a smooth, slowly decaying error.

.. GENERATED FROM PYTHON SOURCE LINES 27-43

.. code-block:: Python


    rng = np.random.default_rng(0)
    n = 33
    e = np.zeros((n, n))
    e[1:-1, 1:-1] = rng.standard_normal((n - 2, n - 2))
    fig1, axes = plt.subplots(1, 3, figsize=(12, 3.8))
    for ax, sweeps in zip(axes, (0, 5, 50), strict=False):
        err = e.copy()
        for _ in range(sweeps):
            err[1:-1, 1:-1] = 0.2 * err[1:-1, 1:-1] + 0.8 * 0.25 * (err[:-2, 1:-1] + err[2:, 1:-1] + err[1:-1, :-2] + err[1:-1, 2:])
        ax.imshow(err.T, origin="lower", cmap="RdBu_r")
        ax.set_title(f"error after {sweeps} Jacobi sweeps")
        ax.set_xticks([])
        ax.set_yticks([])
    fig1.tight_layout()




.. image-sg:: /api/gallery/pde/elliptic/images/sphx_glr_plot_04_multigrid_001.png
   :alt: error after 0 Jacobi sweeps, error after 5 Jacobi sweeps, error after 50 Jacobi sweeps
   :srcset: /api/gallery/pde/elliptic/images/sphx_glr_plot_04_multigrid_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 44-46

Cycle counts stay flat; Jacobi's explode
----------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 46-68

.. code-block:: Python



    def source(X, Y):
        return -2 * np.pi**2 * np.sin(np.pi * X) * np.sin(np.pi * Y)


    for m in (9, 17, 33, 65, 129, 257):
        mg = multigrid_poisson_2d(source, n=m)
        line = f"n = {m:3d} ({(m - 2) ** 2:6d} unknowns): multigrid {mg.solver_result.iterations:2d} V-cycles"
        if m <= 17:
            jac = solve_poisson_2d(source, n=(m, m), method="jacobi", tol=1e-8)
            line += f", Jacobi {jac.solver_result.iterations:4d} sweeps"
        print(line)

    fig2, ax2 = plt.subplots(figsize=(6, 4))
    for m in (17, 65, 257):
        plot_residual_history(multigrid_poisson_2d(source, n=m).solver_result, ax=ax2, label=f"n = {m}")
    ax2.set_xlabel("V-cycle")
    ax2.set_title("Multigrid convergence is independent of grid size")
    fig2.tight_layout()

    plt.show()



.. image-sg:: /api/gallery/pde/elliptic/images/sphx_glr_plot_04_multigrid_002.png
   :alt: Multigrid convergence is independent of grid size
   :srcset: /api/gallery/pde/elliptic/images/sphx_glr_plot_04_multigrid_002.png
   :class: sphx-glr-single-img


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

 .. code-block:: none

    n =   9 (    49 unknowns): multigrid 11 V-cycles, Jacobi  233 sweeps
    n =  17 (   225 unknowns): multigrid 12 V-cycles, Jacobi  950 sweeps
    n =  33 (   961 unknowns): multigrid 12 V-cycles
    n =  65 (  3969 unknowns): multigrid 12 V-cycles
    n = 129 ( 16129 unknowns): multigrid 12 V-cycles
    n = 257 ( 65025 unknowns): multigrid 12 V-cycles





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

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


.. _sphx_glr_download_api_gallery_pde_elliptic_plot_04_multigrid.py:

.. only:: html

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

    .. container:: lite-badge

      .. image:: images/jupyterlite_badge_logo.svg
        :target: ../../../../lite/lab/index.html?path=api/gallery/pde/elliptic/plot_04_multigrid.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

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

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

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

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

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

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


.. only:: html

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

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