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

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

.. _sphx_glr_api_gallery_numerical_analysis_root_finding_plot_01_babylonian_bisection.py:


The Babylonian square root and bisection
========================================

Old Babylonian tablets give :math:`\sqrt{2}` to about six decimal
places. The rule Heron of Alexandria later wrote down averages a guess
with :math:`a` divided by the guess, :math:`x_{n+1} = \tfrac12(x_n + a/x_n)`.
Bisection needs even less: it halves a bracket :math:`[a, b]` on which
:math:`f` changes sign, relying only on the intermediate value theorem.
This script computes :math:`\sqrt{2}` both ways and shows the bracket
shrinking by exactly half per step, against the Babylonian rule's
doubling of correct digits.

.. GENERATED FROM PYTHON SOURCE LINES 16-23

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from mathematicskit.numerical_analysis import Bisection, FixedPointIteration, estimate_convergence_order

    root_exact = np.sqrt(2.0)








.. GENERATED FROM PYTHON SOURCE LINES 24-29

The Babylonian (Heron) iteration
--------------------------------
:math:`g(x) = \tfrac12(x + 2/x)` has :math:`\sqrt{2}` as its fixed point.
Starting from the tablet-friendly guess 1.5, six-place accuracy takes
only three averaging steps.

.. GENERATED FROM PYTHON SOURCE LINES 29-34

.. code-block:: Python


    babylonian = FixedPointIteration(lambda x: 0.5 * (x + 2.0 / x), x0=1.5, tol=1e-15).solve()
    for n, x in enumerate(babylonian.history):
        print(f"Babylonian step {n}: x = {x:.15f}   error = {abs(x - root_exact):.1e}")





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

 .. code-block:: none

    Babylonian step 0: x = 1.500000000000000   error = 8.6e-02
    Babylonian step 1: x = 1.416666666666667   error = 2.5e-03
    Babylonian step 2: x = 1.414215686274510   error = 2.1e-06
    Babylonian step 3: x = 1.414213562374690   error = 1.6e-12
    Babylonian step 4: x = 1.414213562373095   error = 2.2e-16
    Babylonian step 5: x = 1.414213562373095   error = 2.2e-16




.. GENERATED FROM PYTHON SOURCE LINES 35-39

Bisection on :math:`f(x) = x^2 - 2` over :math:`[1, 2]`
-------------------------------------------------------
Each step keeps the half of the bracket where :math:`f` changes sign, so
the bracket width after :math:`n` steps is exactly :math:`2^{-n}`.

.. GENERATED FROM PYTHON SOURCE LINES 39-44

.. code-block:: Python


    bisection = Bisection(lambda x: x**2 - 2.0, 1.0, 2.0, tol=1e-13).solve()
    print(f"bisection: root = {bisection.root:.13f} after {bisection.iterations} halvings")
    print(f"Babylonian order ~ {estimate_convergence_order(babylonian.history[:-1], root_exact):.2f}")





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

 .. code-block:: none

    bisection: root = 1.4142135623730 after 44 halvings
    Babylonian order ~ 1.98




.. GENERATED FROM PYTHON SOURCE LINES 45-47

Bracket width vs. error
-----------------------

.. GENERATED FROM PYTHON SOURCE LINES 47-62

.. code-block:: Python


    fig, ax = plt.subplots(figsize=(7, 4.5))
    steps = np.arange(len(bisection.history))
    bis_err = np.abs(bisection.history - root_exact)
    bab_err = np.abs(babylonian.history - root_exact)
    ax.semilogy(steps, 2.0 ** (-steps.astype(float)), "--", color="gray", label=r"bracket width $2^{-n}$")
    ax.semilogy(steps, np.where(bis_err == 0, np.nan, bis_err), "o-", ms=3, color="darkorange", label="bisection midpoint")
    ax.semilogy(np.arange(len(bab_err)), np.where(bab_err == 0, np.nan, bab_err), "s-", color="seagreen", label="Babylonian average")
    ax.set_xlabel("step $n$")
    ax.set_ylabel(r"$|x_n - \sqrt{2}|$")
    ax.set_title(r"Computing $\sqrt{2}$: halving a bracket vs. Babylonian averaging")
    ax.legend()
    fig.tight_layout()

    plt.show()



.. image-sg:: /api/gallery/numerical_analysis/root_finding/images/sphx_glr_plot_01_babylonian_bisection_001.png
   :alt: Computing $\sqrt{2}$: halving a bracket vs. Babylonian averaging
   :srcset: /api/gallery/numerical_analysis/root_finding/images/sphx_glr_plot_01_babylonian_bisection_001.png
   :class: sphx-glr-single-img






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

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


.. _sphx_glr_download_api_gallery_numerical_analysis_root_finding_plot_01_babylonian_bisection.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_babylonian_bisection.ipynb <plot_01_babylonian_bisection.ipynb>`

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

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

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

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


.. only:: html

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

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