.. 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_03_stormer_verlet.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_03_stormer_verlet.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_03_stormer_verlet.py:


Stormer-Verlet: bounded energy error over long times
====================================================

Stormer (1907) computed charged-particle orbits in the Earth's magnetic
field with the scheme later rediscovered by Verlet (1967) for molecular
dynamics. For :math:`\ddot q = F(q)`:

.. math::

    v_{n+1/2} = v_n + \tfrac{h}{2}F(q_n),\quad
    q_{n+1} = q_n + h\,v_{n+1/2},\quad
    v_{n+1} = v_{n+1/2} + \tfrac{h}{2}F(q_{n+1}).

It is only second order, but it is *symplectic*: it exactly conserves a
slightly perturbed energy. So its energy error stays bounded over
arbitrarily long runs. RK4 is more accurate per step, but its energy
error drifts steadily.

.. GENERATED FROM PYTHON SOURCE LINES 22-42

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from mathematicskit.integrators import leapfrog_integrate, njit, rk4_integrate


    @njit
    def pendulum_force(q, t, params):
        return -np.sin(q)


    @njit
    def pendulum_rhs(state, t, params):
        return np.array([state[1], -np.sin(state[0])])


    def energy(q, v):
        return 0.5 * v**2 - np.cos(q)









.. GENERATED FROM PYTHON SOURCE LINES 43-45

A pendulum released at 2.5 rad, integrated for about 1000 swings
----------------------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 45-64

.. code-block:: Python


    h, n = 0.25, 30_000
    q0, v0 = 2.5, 0.0
    ts, qs, vs = leapfrog_integrate(pendulum_force, np.array([q0]), np.array([v0]), 0.0, h, n, np.zeros(1))
    _, ys = rk4_integrate(pendulum_rhs, np.array([q0, v0]), 0.0, h, n, np.zeros(1))
    e0 = energy(q0, v0)

    fig, ax = plt.subplots(figsize=(8, 4.5))
    ax.plot(ts, energy(ys[:, 0], ys[:, 1]) - e0, lw=0.8, label="RK4 (4th order)")
    ax.plot(ts, energy(qs[:, 0], vs[:, 0]) - e0, lw=0.8, label="Stormer-Verlet (2nd order, symplectic)")
    ax.set_xlabel("t")
    ax.set_ylabel("E(t) - E(0)")
    ax.set_title(f"Pendulum energy error, h = {h}")
    ax.legend()
    fig.tight_layout()

    print(f"max |energy error|: Verlet {np.max(np.abs(energy(qs[:, 0], vs[:, 0]) - e0)):.2e}, RK4 {np.max(np.abs(energy(ys[:, 0], ys[:, 1]) - e0)):.2e}")

    plt.show()



.. image-sg:: /api/gallery/ode_dynamics/integrators/images/sphx_glr_plot_03_stormer_verlet_001.png
   :alt: Pendulum energy error, h = 0.25
   :srcset: /api/gallery/ode_dynamics/integrators/images/sphx_glr_plot_03_stormer_verlet_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none

    max |energy error|: Verlet 1.98e-02, RK4 8.00e-02





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

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


.. _sphx_glr_download_api_gallery_ode_dynamics_integrators_plot_03_stormer_verlet.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_03_stormer_verlet.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

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

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

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

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

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

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


.. only:: html

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

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