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

.. only:: html

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

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

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

.. _sphx_glr_auto_examples_plot_integration_rules.py:


Integration Rules
=================

This example is intended to showcase how integration rules work. There are currently two
supported rules:

- ``"gauss"``, which is uses Gaussian quadrature,
- ``"gauss-lobatto"``, which is uses Gaussian-Lobbato-Legendre nodes for quadrature.

.. GENERATED FROM PYTHON SOURCE LINES 11-16

.. code-block:: Python


    import numpy as np
    from fdg import IntegrationMethod, IntegrationSpecs
    from matplotlib import pyplot as plt








.. GENERATED FROM PYTHON SOURCE LINES 17-25

Different Types of Integration Rules
------------------------------------

There are different types of integration rules. The first one is the Gaussian
quadrature, which is accurate up to the integrals of order :math:`2n - 1` for a rule of
order :math:`n`. The second one is the Gauss-Lobatto quadrature, which also includes the
endpoints of the domain, which makes it useful for some cases, but it reduces the
accuracy of integration to order :math:`2n - 3` for a rule of order :math:`n`.

.. GENERATED FROM PYTHON SOURCE LINES 26-39

.. code-block:: Python


    rule_gauss = IntegrationSpecs(5, IntegrationMethod.GAUSS)
    rule_gl = IntegrationSpecs(5, IntegrationMethod.GAUSS_LOBATTO)

    fig, ax = plt.subplots()
    ax.scatter(rule_gauss.nodes(), rule_gauss.weights(), label="Gauss")
    ax.scatter(rule_gl.nodes(), rule_gl.weights(), label="Gauss-Lobatto")
    ax.grid()
    ax.legend()
    ax.set_ylim(0)
    ax.set(xlabel="$x$", ylabel="$w$")
    plt.show()




.. image-sg:: /auto_examples/images/sphx_glr_plot_integration_rules_001.png
   :alt: plot integration rules
   :srcset: /auto_examples/images/sphx_glr_plot_integration_rules_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 40-46

Integrating Polynomials
-----------------------

Polynomials are exactly integrated up to the degree which is given by the
`IntegrationRule.accuracy` property. The error decreases exponentially
as that order is approached.

.. GENERATED FROM PYTHON SOURCE LINES 47-54

.. code-block:: Python


    rng = np.random.default_rng(seed=0)  # Make rng predictable
    N_POLY = 14
    coeffs = rng.uniform(0, 1, size=N_POLY + 1)
    poly = np.polynomial.Legendre(coeffs, domain=(-1, +1))









.. GENERATED FROM PYTHON SOURCE LINES 55-58

Using Gaussian Quadrature
~~~~~~~~~~~~~~~~~~~~~~~~~


.. GENERATED FROM PYTHON SOURCE LINES 59-73

.. code-block:: Python


    gauss_integral_values: list[float] = list()

    gauss_order = 0
    while True:
        rule = IntegrationSpecs(gauss_order, IntegrationMethod.GAUSS)
        value = np.dot(rule.weights(), poly((rule.nodes() - 1) / 2)) / 2
        gauss_integral_values.append(value)
        if rule.accuracy >= N_POLY + 2:
            # Go one further
            break
        gauss_order += 1









.. GENERATED FROM PYTHON SOURCE LINES 74-77

Using Gaussian-Lobatto Quadrature
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


.. GENERATED FROM PYTHON SOURCE LINES 78-115

.. code-block:: Python


    gauss_lobatto_integral_values: list[float] = list()

    gauss_lobatto_order = 0
    while True:
        rule = IntegrationSpecs(gauss_lobatto_order, IntegrationMethod.GAUSS_LOBATTO)
        value = np.dot(rule.weights(), poly((rule.nodes() - 1) / 2)) / 2
        gauss_lobatto_integral_values.append(value)
        if rule.accuracy >= N_POLY + 2:
            # Go one further
            break
        gauss_lobatto_order += 1


    fig, ax = plt.subplots()

    ax.scatter(
        np.arange(gauss_order),
        np.abs(
            (np.array(gauss_integral_values[:-1]) - gauss_integral_values[-1])
            / gauss_integral_values[-1]
        ),
        label="Gauss",
    )
    ax.scatter(
        np.arange(gauss_lobatto_order),
        np.abs(
            (np.array(gauss_lobatto_integral_values[:-1]) - gauss_lobatto_integral_values[-1])
            / gauss_lobatto_integral_values[-1]
        ),
        label="Gauss-Lobatto",
    )
    ax.set(yscale="log", xlabel="$n$", ylabel="$\\epsilon$")
    ax.grid()
    ax.legend()
    fig.tight_layout()
    plt.show()



.. image-sg:: /auto_examples/images/sphx_glr_plot_integration_rules_002.png
   :alt: plot integration rules
   :srcset: /auto_examples/images/sphx_glr_plot_integration_rules_002.png
   :class: sphx-glr-single-img






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

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


.. _sphx_glr_download_auto_examples_plot_integration_rules.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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