.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/classical/integrators/plot_02_estimate_dt.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. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_classical_integrators_plot_02_estimate_dt.py: Automated timestep selection with estimate_dt =================================================== None of physicskit.classical's integrators are adaptive: each takes a fixed step :math:`\Delta t` for the whole run, so choosing :math:`\Delta t` is normally a manual "integrate, check drift, halve or double, repeat" search -- exactly what :func:`~physicskit.classical.utils.stepsize.estimate_dt` automates, by searching for the largest :math:`\Delta t` whose relative energy drift :math:`|H(t) - H(0)| / |H(0)|` stays below a given tolerance over the full run. This runs it on two very different conservative systems: a quasi-periodic Kepler orbit (:class:`~physicskit.classical.systems.newtonian.KeplerSystem`, semi-major axis :math:`a=1`, eccentricity :math:`e=0.5`, governed by the inverse-square gravitational Hamiltonian), where a short probe already predicts the right :math:`\Delta t`, and a chaotic double pendulum (:class:`~physicskit.classical.systems.lagrangian.DoublePendulum`, governed by its coupled Euler-Lagrange equations), where :func:`estimate_dt`'s own verification stage has to catch and shrink a fast-probe estimate that would otherwise fail once actually run for the full step count. .. GENERATED FROM PYTHON SOURCE LINES 26-33 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from physicskit.classical.systems.lagrangian import DoublePendulum from physicskit.classical.systems.newtonian import KeplerSystem from physicskit.classical.utils.stepsize import estimate_dt .. GENERATED FROM PYTHON SOURCE LINES 34-36 A quasi-periodic system (Kepler orbit) ------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 36-48 .. code-block:: Python print("Quasi-periodic system (Kepler orbit):") est_kepler = estimate_dt( lambda: KeplerSystem.from_orbital_elements(a=1.0, e=0.5), method="yoshida4", tol=1e-6, target_steps=100_000, probe_steps=2000, dt_initial=1e-2, ) print(f" dt = {est_kepler.dt:.4g}, achieved drift = {est_kepler.achieved_drift:.2e}, evaluations = {est_kepler.n_evaluations}") .. rst-class:: sphx-glr-script-out .. code-block:: none Quasi-periodic system (Kepler orbit): dt = 0.01477, achieved drift = 9.92e-07, evaluations = 7 .. GENERATED FROM PYTHON SOURCE LINES 49-55 A chaotic system (double pendulum) --------------------------------------- Found via bracket+bisect on a short probe, THEN verified -- and if necessary shrunk further -- at the full ``target_steps``, which is what makes this safe for a chaotic system where a short probe alone can be misleading. .. GENERATED FROM PYTHON SOURCE LINES 55-67 .. code-block:: Python print("Chaotic system (double pendulum):") est_dp = estimate_dt( lambda: DoublePendulum([2.0, 1.0], [0.5, -0.3]), method="implicit_midpoint", tol=1e-6, target_steps=100_000, probe_steps=2000, dt_initial=1e-3, ) print(f" dt = {est_dp.dt:.4g}, achieved drift = {est_dp.achieved_drift:.2e}, evaluations = {est_dp.n_evaluations}") .. rst-class:: sphx-glr-script-out .. code-block:: none Chaotic system (double pendulum): dt = 3.716e-05, achieved drift = 4.14e-07, evaluations = 11 .. GENERATED FROM PYTHON SOURCE LINES 68-85 Making the search visible: a direct drift-vs-dt sweep ------------------------------------------------------------ The prints above show only the two endpoints of the search -- ``estimate_dt``'s starting guess and its final, verified answer. A direct sweep over dt, all measured on the same short ``probe_steps`` window the fast bracket-and-bisect stage itself uses, draws the underlying drift-vs-dt curve each search is bisecting along -- and makes the difference between the two systems concrete rather than just asserted. For the Kepler orbit, the probe-window curve crosses ``tol`` almost exactly at the ``dt`` ``estimate_dt`` actually returns: the short probe is trustworthy on its own. For the double pendulum, the probe-window curve only crosses ``tol`` at a ``dt`` several times *larger* than ``estimate_dt``'s returned value -- a fast probe alone would have picked a ``dt`` that looks fine over 2000 steps but fails once the chaotic trajectory is actually run out to the full ``target_steps``, which is exactly why the verification-and-shrink stage exists. .. GENERATED FROM PYTHON SOURCE LINES 85-123 .. code-block:: Python def sweep_drift(system_factory, method, dts, n_steps): drifts = np.empty_like(dts) for i, dt in enumerate(dts): system = system_factory() e0 = system.energy() result = system.integrate((0.0, n_steps * dt), dt=dt, method=method) drifts[i] = np.max(np.abs(result.energy - e0)) / abs(e0) return drifts probe_steps = 2000 dts_kepler = np.geomspace(2e-3, 0.2, 12) drift_kepler = sweep_drift(lambda: KeplerSystem.from_orbital_elements(a=1.0, e=0.5), "yoshida4", dts_kepler, probe_steps) dts_dp = np.geomspace(1e-5, 3e-3, 12) drift_dp = sweep_drift(lambda: DoublePendulum([2.0, 1.0], [0.5, -0.3]), "implicit_midpoint", dts_dp, probe_steps) fig, axes = plt.subplots(1, 2, figsize=(10, 4.3)) axes[0].loglog(dts_kepler, drift_kepler, "o-", color="steelblue") axes[0].axhline(1e-6, color="0.6", ls="--", lw=1, label="tol") axes[0].axvline(est_kepler.dt, color="firebrick", ls=":", lw=1.2, label=f"estimate_dt: {est_kepler.dt:.3g}") axes[0].set_xlabel("dt") axes[0].set_ylabel("relative energy drift (probe window)") axes[0].set_title("Kepler orbit: probe crossing = estimate_dt's answer", fontsize=10.5) axes[0].legend(fontsize=8) axes[1].loglog(dts_dp, drift_dp, "o-", color="firebrick") axes[1].axhline(1e-6, color="0.6", ls="--", lw=1, label="tol") axes[1].axvline(est_dp.dt, color="steelblue", ls=":", lw=1.2, label=f"estimate_dt: {est_dp.dt:.3g}") axes[1].set_xlabel("dt") axes[1].set_ylabel("relative energy drift (probe window)") axes[1].set_title("Double pendulum: probe crossing overshoots the answer", fontsize=10.5) axes[1].legend(fontsize=8) fig.tight_layout() plt.show() .. image-sg:: /api/gallery/classical/integrators/images/sphx_glr_plot_02_estimate_dt_001.png :alt: Kepler orbit: probe crossing = estimate_dt's answer, Double pendulum: probe crossing overshoots the answer :srcset: /api/gallery/classical/integrators/images/sphx_glr_plot_02_estimate_dt_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 5.354 seconds) .. _sphx_glr_download_api_gallery_classical_integrators_plot_02_estimate_dt.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_02_estimate_dt.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_02_estimate_dt.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_02_estimate_dt.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_