.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/calculus/quadrature/plot_01_gauss_legendre.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_calculus_quadrature_plot_01_gauss_legendre.py>`
        to download the full example code.

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

.. _sphx_glr_api_gallery_calculus_quadrature_plot_01_gauss_legendre.py:


Gaussian quadrature: optimal nodes are exact to degree 2n-1
=============================================================

Gauss let a quadrature rule choose its sample points as well as its
weights. The optimal :math:`n` nodes turn out to be the roots of the
Legendre polynomial :math:`P_n`, and the resulting rule integrates every
polynomial of degree up to :math:`2n-1` exactly -- about twice the degree
an equally spaced rule with the same number of points manages. This
example shows the nodes and weights from
:func:`~mathematicskit.calculus.legendre_nodes_and_weights`, checks the
degree of exactness of
:class:`~mathematicskit.calculus.GaussianQuadrature`, and compares its
convergence with the equally spaced trapezoidal rule.

.. GENERATED FROM PYTHON SOURCE LINES 18-24

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np
    from numpy.polynomial import legendre

    from mathematicskit.calculus import GaussianQuadrature, TrapezoidalRule, legendre_nodes_and_weights








.. GENERATED FROM PYTHON SOURCE LINES 25-27

Nodes are the roots of the Legendre polynomial
----------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 27-45

.. code-block:: Python


    n = 5
    nodes, weights = legendre_nodes_and_weights(n)
    grid = np.linspace(-1, 1, 400)
    p_n = legendre.Legendre.basis(n)

    fig, ax = plt.subplots(figsize=(6, 4))
    ax.plot(grid, p_n(grid), "k", label=f"Legendre $P_{n}$")
    ax.axhline(0, color="gray", lw=0.5)
    ax.stem(nodes, weights, linefmt="C0-", markerfmt="C0o", basefmt=" ", label="Gauss nodes (height = weight)")
    ax.set_title(f"{n}-point Gauss-Legendre rule on [-1, 1]")
    ax.legend()
    fig.tight_layout()

    print("nodes:  ", np.round(nodes, 6))
    print("weights:", np.round(weights, 6))
    print("max |P_n(node)|:", float(np.max(np.abs(p_n(nodes)))))




.. image-sg:: /api/gallery/calculus/quadrature/images/sphx_glr_plot_01_gauss_legendre_001.png
   :alt: 5-point Gauss-Legendre rule on [-1, 1]
   :srcset: /api/gallery/calculus/quadrature/images/sphx_glr_plot_01_gauss_legendre_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none

    nodes:   [-0.90618  -0.538469  0.        0.538469  0.90618 ]
    weights: [0.236927 0.478629 0.568889 0.478629 0.236927]
    max |P_n(node)|: 1.1102230246251565e-16




.. GENERATED FROM PYTHON SOURCE LINES 46-51

Exact for every polynomial of degree up to 2n - 1
-------------------------------------------------

With :math:`n = 5` nodes, :math:`x^k` over :math:`[0, 1]` is integrated
to round-off for :math:`k \le 9` and first fails at :math:`k = 10`.

.. GENERATED FROM PYTHON SOURCE LINES 51-57

.. code-block:: Python


    rule = GaussianQuadrature(n=n)
    for k in range(2 * n + 2):
        err = abs(rule.integrate(lambda x, k=k: x**k, 0.0, 1.0).value - 1.0 / (k + 1))
        print(f"degree {k:2d}: error = {err:.2e}{'  <- exactness lost' if k == 2 * n else ''}")





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

 .. code-block:: none

    degree  0: error = 0.00e+00
    degree  1: error = 0.00e+00
    degree  2: error = 5.55e-17
    degree  3: error = 1.11e-16
    degree  4: error = 8.33e-17
    degree  5: error = 8.33e-17
    degree  6: error = 1.11e-16
    degree  7: error = 8.33e-17
    degree  8: error = 9.71e-17
    degree  9: error = 8.33e-17
    degree 10: error = 1.43e-06  <- exactness lost
    degree 11: error = 7.87e-06




.. GENERATED FROM PYTHON SOURCE LINES 58-60

Far fewer evaluations than an equally spaced rule
-------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 60-76

.. code-block:: Python


    f, a, b, exact = np.exp, 0.0, 1.0, np.e - 1.0
    ns = np.arange(1, 11)
    err_gauss = [abs(GaussianQuadrature(n=int(k)).integrate(f, a, b).value - exact) for k in ns]
    err_trap = [abs(TrapezoidalRule(n=int(k) - 1).integrate(f, a, b).value - exact) if k > 1 else np.nan for k in ns]

    fig, ax = plt.subplots(figsize=(6, 4.5))
    ax.semilogy(ns, np.maximum(err_gauss, 1e-17), "o-", label="Gauss-Legendre")
    ax.semilogy(ns, err_trap, "s-", label="trapezoidal (equally spaced)")
    ax.set_xlabel("function evaluations")
    ax.set_ylabel("|error|")
    ax.set_title("Integral of exp over [0, 1]")
    ax.legend()
    fig.tight_layout()

    plt.show()



.. image-sg:: /api/gallery/calculus/quadrature/images/sphx_glr_plot_01_gauss_legendre_002.png
   :alt: Integral of exp over [0, 1]
   :srcset: /api/gallery/calculus/quadrature/images/sphx_glr_plot_01_gauss_legendre_002.png
   :class: sphx-glr-single-img






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

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


.. _sphx_glr_download_api_gallery_calculus_quadrature_plot_01_gauss_legendre.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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