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

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

.. _sphx_glr_api_gallery_special_functions_elliptic_plot_01_agm_elliptic_integrals.py:


Gauss's arithmetic-geometric mean and the elliptic integral K
===================================================================

Iterates the AGM by hand to show its quadratic convergence, then checks
Gauss's 1799 identity K(m) = pi / (2 AGM(1, sqrt(1 - m))) and uses E(m)
for the perimeter of an ellipse.

.. GENERATED FROM PYTHON SOURCE LINES 11-18

.. code-block:: Python

    import math

    import matplotlib.pyplot as plt
    import numpy as np

    from mathematicskit.special_functions import arithmetic_geometric_mean, complete_elliptic_integral_first_kind, complete_elliptic_integral_second_kind








.. GENERATED FROM PYTHON SOURCE LINES 19-21

Quadratic convergence: the number of correct digits doubles each step
----------------------------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 21-29

.. code-block:: Python


    a, b = 1.0, math.sqrt(2.0)
    limit = arithmetic_geometric_mean(a, b)
    for step in range(5):
        print(f"step {step}: a = {a:.16f}, b = {b:.16f}, |a - AGM| = {abs(a - limit):.1e}")
        a, b = (a + b) / 2, math.sqrt(a * b)
    print(f"Gauss's constant 1/AGM(1, sqrt 2) = {1 / limit:.12f}")





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

 .. code-block:: none

    step 0: a = 1.0000000000000000, b = 1.4142135623730951, |a - AGM| = 2.0e-01
    step 1: a = 1.2071067811865475, b = 1.1892071150027210, |a - AGM| = 9.0e-03
    step 2: a = 1.1981569480946344, b = 1.1981235214931201, |a - AGM| = 1.7e-05
    step 3: a = 1.1981402347938772, b = 1.1981402346773073, |a - AGM| = 5.8e-11
    step 4: a = 1.1981402347355923, b = 1.1981402347355923, |a - AGM| = 2.2e-16
    Gauss's constant 1/AGM(1, sqrt 2) = 0.834626841674




.. GENERATED FROM PYTHON SOURCE LINES 30-32

K(m) from the AGM
-----------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 32-45

.. code-block:: Python


    m = np.linspace(0.0, 0.999, 200)
    k_agm = np.pi / (2 * arithmetic_geometric_mean(1.0, np.sqrt(1 - m)))
    print(f"\nmax |K(m) - pi/(2 AGM)| = {np.max(np.abs(complete_elliptic_integral_first_kind(m) - k_agm)):.1e}")

    fig, ax = plt.subplots()
    ax.plot(m, complete_elliptic_integral_first_kind(m), label="K(m)")
    ax.plot(m, complete_elliptic_integral_second_kind(m), label="E(m)")
    ax.axhline(np.pi / 2, color="gray", ls=":", lw=1)
    ax.set_xlabel("parameter m")
    ax.set_title("Complete elliptic integrals")
    ax.legend()




.. image-sg:: /api/gallery/special_functions/elliptic/images/sphx_glr_plot_01_agm_elliptic_integrals_001.png
   :alt: Complete elliptic integrals
   :srcset: /api/gallery/special_functions/elliptic/images/sphx_glr_plot_01_agm_elliptic_integrals_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none


    max |K(m) - pi/(2 AGM)| = 8.9e-16

    <matplotlib.legend.Legend object at 0x11a99dd30>



.. GENERATED FROM PYTHON SOURCE LINES 46-48

The perimeter of an ellipse is 4 a E(e^2)
-----------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 48-55

.. code-block:: Python


    semi_major, semi_minor = 2.0, 1.0
    ecc2 = 1 - (semi_minor / semi_major) ** 2
    perimeter = 4 * semi_major * complete_elliptic_integral_second_kind(ecc2)
    t = np.linspace(0, 2 * np.pi, 200001)
    polyline = np.sum(np.hypot(np.diff(semi_major * np.cos(t)), np.diff(semi_minor * np.sin(t))))
    print(f"\nellipse 2 x 1: 4aE(e^2) = {perimeter:.10f}, polyline length = {polyline:.10f}")




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

 .. code-block:: none


    ellipse 2 x 1: 4aE(e^2) = 9.6884482205, polyline length = 9.6884482201





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

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


.. _sphx_glr_download_api_gallery_special_functions_elliptic_plot_01_agm_elliptic_integrals.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_agm_elliptic_integrals.ipynb <plot_01_agm_elliptic_integrals.ipynb>`

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

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

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

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


.. only:: html

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

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