.. 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_05_yoshida_composition.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` 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_05_yoshida_composition.py: Yoshida's symplectic composition: fourth order from leapfrog ============================================================ Yoshida (1990) raised the order of a symmetric method by *composing* it with itself. Three leapfrog substeps of sizes :math:`w_1 h, w_0 h, w_1 h` with .. math:: w_1 = \frac{1}{2 - 2^{1/3}},\qquad w_0 = -\frac{2^{1/3}}{2 - 2^{1/3}} cancel the :math:`h^3` error term. The result is fourth order and, being built from symplectic steps, still symplectic. The middle substep runs *backwards* in time (:math:`w_0 < 0`). .. GENERATED FROM PYTHON SOURCE LINES 19-39 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from mathematicskit.integrators import leapfrog_integrate, njit, yoshida4_integrate @njit def spring(q, t, params): return -q t_end = 10.0 ns = 2 ** np.arange(5, 12) hs = t_end / ns errors = {"leapfrog": [], "Yoshida 4th order": []} for h, n in zip(hs, ns, strict=False): for label, method in (("leapfrog", leapfrog_integrate), ("Yoshida 4th order", yoshida4_integrate)): _, qs, _ = method(spring, np.array([1.0]), np.array([0.0]), 0.0, h, n, np.zeros(1)) errors[label].append(abs(qs[-1, 0] - np.cos(t_end))) .. GENERATED FROM PYTHON SOURCE LINES 40-42 Error vs step size ------------------ .. GENERATED FROM PYTHON SOURCE LINES 42-56 .. code-block:: Python fig, ax = plt.subplots(figsize=(7, 4.5)) for (label, err), cost in zip(errors.items(), (1, 3), strict=False): ax.loglog(hs, err, "o-", label=f"{label} ({cost} force evaluation(s) per step)") print(f"{label}: error ratio for halved h = {err[-2] / err[-1]:.1f}") ax.loglog(hs, 0.05 * hs**2, "k--", lw=0.8, label=r"$O(h^2)$") ax.loglog(hs, 0.01 * hs**4, "k:", lw=0.8, label=r"$O(h^4)$") ax.set_xlabel("h") ax.set_ylabel("|q(10) - cos 10|") ax.set_title("Yoshida composition: slope 4 instead of 2") ax.legend() fig.tight_layout() plt.show() .. image-sg:: /api/gallery/ode_dynamics/integrators/images/sphx_glr_plot_05_yoshida_composition_001.png :alt: Yoshida composition: slope 4 instead of 2 :srcset: /api/gallery/ode_dynamics/integrators/images/sphx_glr_plot_05_yoshida_composition_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none leapfrog: error ratio for halved h = 4.0 Yoshida 4th order: error ratio for halved h = 16.0 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 2.103 seconds) .. _sphx_glr_download_api_gallery_ode_dynamics_integrators_plot_05_yoshida_composition.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_05_yoshida_composition.ipynb :alt: Launch JupyterLite :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_05_yoshida_composition.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_05_yoshida_composition.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_05_yoshida_composition.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_