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

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

.. _sphx_glr_api_gallery_chaos_chaos_metrics_plot_bifurcation_diagram.py:


Bifurcation Diagrams for Continuous Flows: the Duffing Oscillator
======================================================================

Bifurcation diagrams aren't only for discrete maps -- they're just as useful
for continuous, periodically-forced flows, via a *stroboscopic* Poincare
section: sampling the state once per forcing period turns the continuous
flow into an effective discrete map. The flow swept here is the forced,
damped Duffing oscillator,

.. math::

    \ddot{x} + \delta \dot{x} + \alpha x + \beta x^3 = \gamma \cos(\omega t),

with fixed damping :math:`\delta=0.3`, linear stiffness :math:`\alpha=-1`,
cubic stiffness :math:`\beta=1`, and forcing frequency :math:`\omega=1.2`
(so the double-well potential and forcing frequency are held fixed while
only the forcing amplitude :math:`\gamma` is swept). This example sweeps
:math:`\gamma` and, for each value, plots the surviving stroboscopic
samples, using
:func:`physicskit.chaos.visualizers.bifurcation.stroboscopic_bifurcation_sampler` --
the same :func:`~physicskit.chaos.visualizers.bifurcation.plot_bifurcation_diagram`
used for the :doc:`Logistic Map </api/gallery/chaos/maps/plot_logistic_map>`
works here completely unchanged. A companion plot below sweeps the same
:math:`\gamma` range with
:func:`physicskit.chaos.visualizers.divergence.trajectory_divergence` and
:func:`physicskit.chaos.utils.metrics.lyapunov_exponent_from_divergence` to
estimate the largest Lyapunov exponent at each amplitude, turning the
bifurcation diagram's purely qualitative "band looks chaotic" into the
quantitative :math:`\lambda_{\max} > 0` criterion, on the same :math:`\gamma`
axis.

.. GENERATED FROM PYTHON SOURCE LINES 33-45

.. code-block:: Python


    import matplotlib.pyplot as plt
    import numpy as np

    from physicskit.chaos.systems.continuous import Duffing
    from physicskit.chaos.utils.metrics import lyapunov_exponent_from_divergence
    from physicskit.chaos.visualizers.bifurcation import (
        plot_bifurcation_diagram,
        stroboscopic_bifurcation_sampler,
    )
    from physicskit.chaos.visualizers.divergence import trajectory_divergence








.. GENERATED FROM PYTHON SOURCE LINES 46-50

Build the stroboscopic sampler
----------------------------------
The forcing period is ``2*pi / omega``; sampling the state once per forcing
period is what turns the continuous flow into an effective discrete map.

.. GENERATED FROM PYTHON SOURCE LINES 50-65

.. code-block:: Python

    omega = 1.2
    period = 2.0 * np.pi / omega

    sampler = stroboscopic_bifurcation_sampler(
        lambda gamma: Duffing(delta=0.3, alpha=-1.0, beta=1.0, gamma=gamma, omega=omega),
        state0=np.array([0.1, 0.0]),
        sample_period=period,
        n_transient_periods=60,
        n_keep_periods=20,
        component=0,
        dt=0.15,
        rtol=1e-6,
        atol=1e-8,
    )








.. GENERATED FROM PYTHON SOURCE LINES 66-77

Sweep the forcing amplitude
--------------------------------
As ``gamma`` increases, watch the single-point regular response give way
to period-doubled cycles, and eventually a chaotic band -- the same
period-doubling route to chaos as the Logistic Map, now in a physically
driven mechanical oscillator. (Looser tolerances than the sampler's
defaults keep this parameter sweep fast -- it calls the integrator once per
`gamma` value below.) Each `gamma` value is fully independent of every
other, so this sweep is also embarrassingly parallel: ``n_jobs`` runs it
across a thread pool (real speedup here, since `solve_ivp` releases the
GIL while it integrates), and ``show_progress`` displays a `tqdm` bar.

.. GENERATED FROM PYTHON SOURCE LINES 77-86

.. code-block:: Python

    gamma_values = np.linspace(0.2, 0.5, 150)
    fig, ax = plt.subplots(figsize=(9, 5))
    plot_bifurcation_diagram(gamma_values, sampler, n_jobs=4, show_progress=True, marker=".", markersize=1.0, ax=ax)
    ax.set_xlabel(r"forcing amplitude $\gamma$")
    ax.set_ylabel("x (stroboscopic samples)")
    ax.set_title("Duffing oscillator: bifurcation diagram vs. forcing amplitude")

    plt.show()




.. image-sg:: /api/gallery/chaos/chaos_metrics/images/sphx_glr_plot_bifurcation_diagram_001.png
   :alt: Duffing oscillator: bifurcation diagram vs. forcing amplitude
   :srcset: /api/gallery/chaos/chaos_metrics/images/sphx_glr_plot_bifurcation_diagram_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none


    sweeping:   0%|          | 0/150 [00:00<?, ?it/s]
    sweeping:   1%|          | 1/150 [00:00<01:18,  1.90it/s]
    sweeping:   3%|▎         | 5/150 [00:00<00:19,  7.41it/s]
    sweeping:   4%|▍         | 6/150 [00:00<00:19,  7.39it/s]
    sweeping:   5%|▌         | 8/150 [00:01<00:17,  7.96it/s]
    sweeping:   7%|▋         | 10/150 [00:01<00:17,  8.21it/s]
    sweeping:   7%|▋         | 11/150 [00:01<00:17,  8.13it/s]
    sweeping:   8%|▊         | 12/150 [00:01<00:17,  7.78it/s]
    sweeping:   9%|▉         | 14/150 [00:01<00:15,  8.87it/s]
    sweeping:  10%|█         | 15/150 [00:02<00:17,  7.87it/s]
    sweeping:  11%|█▏        | 17/150 [00:02<00:14,  8.94it/s]
    sweeping:  13%|█▎        | 19/150 [00:02<00:13,  9.56it/s]
    sweeping:  13%|█▎        | 20/150 [00:02<00:14,  9.03it/s]
    sweeping:  14%|█▍        | 21/150 [00:02<00:16,  7.67it/s]
    sweeping:  16%|█▌        | 24/150 [00:03<00:15,  8.32it/s]
    sweeping:  17%|█▋        | 25/150 [00:03<00:16,  7.46it/s]
    sweeping:  19%|█▊        | 28/150 [00:03<00:14,  8.19it/s]
    sweeping:  21%|██        | 31/150 [00:03<00:11, 10.40it/s]
    sweeping:  22%|██▏       | 33/150 [00:04<00:13,  8.63it/s]
    sweeping:  23%|██▎       | 35/150 [00:04<00:12,  9.00it/s]
    sweeping:  25%|██▍       | 37/150 [00:04<00:13,  8.16it/s]
    sweeping:  26%|██▌       | 39/150 [00:04<00:13,  8.12it/s]
    sweeping:  27%|██▋       | 41/150 [00:04<00:11,  9.29it/s]
    sweeping:  29%|██▊       | 43/150 [00:05<00:09, 10.75it/s]
    sweeping:  30%|███       | 45/150 [00:05<00:12,  8.57it/s]
    sweeping:  31%|███▏      | 47/150 [00:05<00:10, 10.26it/s]
    sweeping:  33%|███▎      | 49/150 [00:05<00:12,  8.29it/s]
    sweeping:  34%|███▍      | 51/150 [00:05<00:10,  9.49it/s]
    sweeping:  35%|███▌      | 53/150 [00:06<00:10,  8.94it/s]
    sweeping:  37%|███▋      | 55/150 [00:06<00:10,  8.65it/s]
    sweeping:  38%|███▊      | 57/150 [00:06<00:10,  8.86it/s]
    sweeping:  39%|███▊      | 58/150 [00:06<00:10,  8.41it/s]
    sweeping:  39%|███▉      | 59/150 [00:06<00:11,  8.11it/s]
    sweeping:  41%|████      | 61/150 [00:07<00:09,  9.68it/s]
    sweeping:  42%|████▏     | 63/150 [00:07<00:09,  9.43it/s]
    sweeping:  43%|████▎     | 64/150 [00:07<00:09,  9.12it/s]
    sweeping:  43%|████▎     | 65/150 [00:07<00:09,  9.21it/s]
    sweeping:  45%|████▍     | 67/150 [00:07<00:09,  8.64it/s]
    sweeping:  45%|████▌     | 68/150 [00:07<00:09,  8.33it/s]
    sweeping:  47%|████▋     | 71/150 [00:08<00:09,  8.30it/s]
    sweeping:  48%|████▊     | 72/150 [00:08<00:09,  8.18it/s]
    sweeping:  50%|█████     | 75/150 [00:08<00:08,  8.71it/s]
    sweeping:  51%|█████▏    | 77/150 [00:08<00:07, 10.34it/s]
    sweeping:  53%|█████▎    | 79/150 [00:09<00:09,  7.71it/s]
    sweeping:  55%|█████▍    | 82/150 [00:09<00:06,  9.87it/s]
    sweeping:  56%|█████▌    | 84/150 [00:09<00:07,  8.65it/s]
    sweeping:  57%|█████▋    | 86/150 [00:10<00:07,  8.84it/s]
    sweeping:  59%|█████▊    | 88/150 [00:10<00:07,  8.31it/s]
    sweeping:  61%|██████    | 91/150 [00:10<00:06,  9.41it/s]
    sweeping:  62%|██████▏   | 93/150 [00:10<00:06,  8.68it/s]
    sweeping:  63%|██████▎   | 95/150 [00:10<00:05,  9.57it/s]
    sweeping:  65%|██████▍   | 97/150 [00:11<00:05,  9.24it/s]
    sweeping:  65%|██████▌   | 98/150 [00:11<00:05,  9.21it/s]
    sweeping:  66%|██████▌   | 99/150 [00:11<00:05,  8.52it/s]
    sweeping:  67%|██████▋   | 101/150 [00:11<00:05,  9.48it/s]
    sweeping:  68%|██████▊   | 102/150 [00:11<00:05,  8.02it/s]
    sweeping:  69%|██████▉   | 104/150 [00:11<00:04,  9.68it/s]
    sweeping:  71%|███████   | 106/150 [00:12<00:04,  9.35it/s]
    sweeping:  71%|███████▏  | 107/150 [00:12<00:05,  8.27it/s]
    sweeping:  73%|███████▎  | 109/150 [00:12<00:04,  9.27it/s]
    sweeping:  74%|███████▍  | 111/150 [00:12<00:05,  7.80it/s]
    sweeping:  75%|███████▌  | 113/150 [00:12<00:03,  9.51it/s]
    sweeping:  77%|███████▋  | 115/150 [00:13<00:04,  7.63it/s]
    sweeping:  79%|███████▉  | 119/150 [00:13<00:03,  9.20it/s]
    sweeping:  81%|████████  | 121/150 [00:13<00:02,  9.85it/s]
    sweeping:  82%|████████▏ | 123/150 [00:14<00:02,  9.50it/s]
    sweeping:  83%|████████▎ | 125/150 [00:14<00:02,  9.28it/s]
    sweeping:  84%|████████▍ | 126/150 [00:14<00:02,  9.03it/s]
    sweeping:  85%|████████▌ | 128/150 [00:14<00:02,  7.46it/s]
    sweeping:  87%|████████▋ | 130/150 [00:14<00:02,  8.67it/s]
    sweeping:  88%|████████▊ | 132/150 [00:15<00:02,  8.64it/s]
    sweeping:  89%|████████▊ | 133/150 [00:15<00:01,  8.52it/s]
    sweeping:  90%|█████████ | 135/150 [00:15<00:01, 10.15it/s]
    sweeping:  91%|█████████▏| 137/150 [00:15<00:01, 10.37it/s]
    sweeping:  93%|█████████▎| 139/150 [00:15<00:01,  8.62it/s]
    sweeping:  94%|█████████▍| 141/150 [00:16<00:00,  9.93it/s]
    sweeping:  95%|█████████▌| 143/150 [00:16<00:00,  7.60it/s]
    sweeping:  97%|█████████▋| 146/150 [00:16<00:00,  7.69it/s]
    sweeping: 100%|██████████| 150/150 [00:16<00:00,  8.84it/s]




.. GENERATED FROM PYTHON SOURCE LINES 87-100

Quantifying it: the largest Lyapunov exponent, on the same axis
---------------------------------------------------------------------
The bifurcation diagram above only shows *where* the response looks
irregular; it doesn't say whether that irregularity is genuine exponential
sensitivity or "just" a very long periodic cycle. Sweeping the same
`gamma_values` (on a coarser grid, since each point here integrates two
whole trajectories rather than sampling one) with
:func:`~physicskit.chaos.visualizers.divergence.trajectory_divergence` and
fitting the growth rate of two initially nearby trajectories'
separation gives the largest Lyapunov exponent :math:`\lambda_{\max}`
directly: it should track the bifurcation diagram closely, crossing zero
right around where the single point (or handful of stroboscopic points)
gives way to the dense chaotic band above.

.. GENERATED FROM PYTHON SOURCE LINES 100-118

.. code-block:: Python

    gamma_values_lyap = np.linspace(0.2, 0.5, 50)
    lyapunov_exponents = np.empty_like(gamma_values_lyap)
    for i, gamma in enumerate(gamma_values_lyap):
        lyap_system = Duffing(delta=0.3, alpha=-1.0, beta=1.0, gamma=float(gamma), omega=omega)
        t_div, delta_t = trajectory_divergence(lyap_system, state0=np.array([0.1, 0.0]), delta_0=1e-8, t_max=150.0, n_points=400, seed=0)
        n_fit = int(0.5 * t_div.size)
        lyapunov_exponents[i] = lyapunov_exponent_from_divergence(t_div[:n_fit], delta_t[:n_fit] / 1e-8)

    fig_lyap, ax_lyap = plt.subplots(figsize=(9, 4))
    ax_lyap.axhline(0.0, color="black", lw=0.8)
    ax_lyap.plot(gamma_values_lyap, lyapunov_exponents, "o-", color="crimson", markersize=3)
    ax_lyap.set_xlim(gamma_values[0], gamma_values[-1])
    ax_lyap.set_xlabel(r"forcing amplitude $\gamma$")
    ax_lyap.set_ylabel(r"largest Lyapunov exponent $\lambda_{max}$")
    ax_lyap.set_title(r"$\lambda_{max} > 0$ lines up with the bifurcation diagram's chaotic band")
    fig_lyap.tight_layout()

    plt.show()



.. image-sg:: /api/gallery/chaos/chaos_metrics/images/sphx_glr_plot_bifurcation_diagram_002.png
   :alt: $\lambda_{max} > 0$ lines up with the bifurcation diagram's chaotic band
   :srcset: /api/gallery/chaos/chaos_metrics/images/sphx_glr_plot_bifurcation_diagram_002.png
   :class: sphx-glr-single-img






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

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


.. _sphx_glr_download_api_gallery_chaos_chaos_metrics_plot_bifurcation_diagram.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

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

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

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

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

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

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


.. only:: html

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

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