.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/ode_dynamics/integrators/plot_04_dormand_prince.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_ode_dynamics_integrators_plot_04_dormand_prince.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_ode_dynamics_integrators_plot_04_dormand_prince.py:


Dormand-Prince: embedded pairs and adaptive step size
=====================================================

An *embedded* Runge-Kutta pair computes a 5th- and a 4th-order solution
from the same seven stages. Their difference estimates the local error
at almost no extra cost. Dormand and Prince (1980) tuned the pair so that
the 5th-order result, the one actually kept, has a very small error
constant. A controller then shrinks :math:`h` where the error estimate is
large and grows it where it is small. On a highly eccentric Kepler orbit
the step collapses at each close pass (perihelion) and stretches out far
from the Sun (aphelion).

.. GENERATED FROM PYTHON SOURCE LINES 16-32

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from mathematicskit.integrators import dopri5_integrate, njit


    @njit
    def kepler(state, t, params):
        r3 = (state[0] ** 2 + state[1] ** 2) ** 1.5
        return np.array([state[2], state[3], -state[0] / r3, -state[1] / r3])


    e = 0.9
    state0 = np.array([1.0 - e, 0.0, 0.0, np.sqrt((1 + e) / (1 - e))])  # perihelion, semi-major axis 1
    period = 2 * np.pi








.. GENERATED FROM PYTHON SOURCE LINES 33-35

Steps crowd around perihelion
-----------------------------

.. GENERATED FROM PYTHON SOURCE LINES 35-52

.. code-block:: Python


    ts, ys = dopri5_integrate(kepler, state0, 0.0, 3 * period, 1e-3, np.zeros(1), rtol=1e-8, atol=1e-10)
    r = np.hypot(ys[:, 0], ys[:, 1])

    fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(11, 4.5))
    ax0.plot(ys[:, 0], ys[:, 1], ".", ms=2)
    ax0.plot(0, 0, "*", color="orange", ms=12)
    ax0.set_aspect("equal")
    ax0.set_title(f"Accepted steps on an e = {e} orbit ({len(ts) - 1} steps)")

    ax1.semilogy(ts[1:], np.diff(ts), label="step size h")
    ax1.semilogy(ts, r, label="distance r")
    ax1.set_xlabel("t")
    ax1.set_title("h tracks the orbit's time scale")
    ax1.legend()
    fig.tight_layout()




.. image-sg:: /api/gallery/ode_dynamics/integrators/images/sphx_glr_plot_04_dormand_prince_001.png
   :alt: Accepted steps on an e = 0.9 orbit (435 steps), h tracks the orbit's time scale
   :srcset: /api/gallery/ode_dynamics/integrators/images/sphx_glr_plot_04_dormand_prince_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 53-55

Tolerance sets the cost
-----------------------

.. GENERATED FROM PYTHON SOURCE LINES 55-61

.. code-block:: Python


    for rtol in (1e-4, 1e-6, 1e-8, 1e-10):
        ts_, ys_ = dopri5_integrate(kepler, state0, 0.0, period, 1e-3, np.zeros(1), rtol=rtol, atol=rtol * 1e-2)
        print(f"rtol={rtol:.0e}: {len(ts_) - 1:5d} steps, error after one period = {np.linalg.norm(ys_[-1] - state0):.2e}")

    plt.show()




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

 .. code-block:: none

    rtol=1e-04:    28 steps, error after one period = 1.41e+00
    rtol=1e-06:    61 steps, error after one period = 1.97e-02
    rtol=1e-08:   145 steps, error after one period = 1.22e-04
    rtol=1e-10:   364 steps, error after one period = 7.02e-07





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

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


.. _sphx_glr_download_api_gallery_ode_dynamics_integrators_plot_04_dormand_prince.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/ode_dynamics/integrators/plot_04_dormand_prince.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

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

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

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

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

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

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


.. only:: html

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

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