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

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

.. _sphx_glr_api_gallery_statphys_stochastic_processes_plot_langevin_fluctuation_dissipation.py:


The Langevin equation: friction and noise must balance
=========================================================

Langevin's 1908 equation for a Brownian particle,

.. math::

    m\,dv = (F(x) - \gamma v)\,dt + \sqrt{2\gamma k_BT}\,dW,

contains a drag and a random kick with the same :math:`\gamma` in both.
This is not a modeling convenience. It is the fluctuation-dissipation
balance: the drag drains energy at a rate proportional to :math:`\gamma`,
the noise pumps it in at a rate proportional to its variance, and they
cancel only at :math:`\langle v^2\rangle = k_BT/m`.

This example integrates the equation for particles in a harmonic trap
with :class:`~physicskit.statphys.LangevinDynamics`, which uses the BAOAB
splitting (:func:`~physicskit.integrators.baoab_integrate`) so that even
weak friction does not let the time step heat the particles. The same temperature
comes out for frictions spanning two decades, and every quadratic degree
of freedom carries :math:`k_BT/2` (equipartition). It then breaks the
balance by integrating the equation with noise that is too weak, using the
shared :func:`~physicskit.integrators.euler_maruyama_integrate` directly,
and the particles settle at the wrong temperature.

.. GENERATED FROM PYTHON SOURCE LINES 27-37

.. code-block:: Python


    import matplotlib.pyplot as plt
    import numpy as np
    from numba import njit

    from physicskit.integrators import euler_maruyama_integrate, wiener_increments
    from physicskit.statphys import LangevinDynamics

    m, k, kT = 1.0, 4.0, 0.5








.. GENERATED FROM PYTHON SOURCE LINES 38-45

One temperature, any friction
-----------------------------
Particles start at rest at the trap's edge. Light friction lets them ring
for many periods and heavy friction makes them creep, yet all settle at
:math:`m\langle v^2\rangle = k\langle x^2\rangle = k_BT`. We average
over :math:`t > 60`, many relaxation times :math:`m/\gamma` even for the
lightest friction.

.. GENERATED FROM PYTHON SOURCE LINES 45-56

.. code-block:: Python


    results = {}
    equipartition = {}
    for gamma in (0.1, 1.0, 10.0):
        ld = LangevinDynamics(n_particles=4000, mass=m, gamma=gamma, kT=kT, stiffness=k, seed=int(10 * gamma))
        t, x, v = ld.run(t_max=80.0, dt=0.002, n_frames=400, x0=np.full(4000, 1.0))
        results[gamma] = (t, x[:, :, 0], v[:, :, 0])
        late = t > 60
        equipartition[gamma] = (m * np.mean(v[late] ** 2), k * np.mean(x[late] ** 2))
        print(f"gamma = {gamma:5.1f}: m<v^2> = {m * np.mean(v[late] ** 2):.3f}, k<x^2> = {k * np.mean(x[late] ** 2):.3f}   (kT = {kT})")





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

 .. code-block:: none

    gamma =   0.1: m<v^2> = 0.499, k<x^2> = 0.500   (kT = 0.5)
    gamma =   1.0: m<v^2> = 0.498, k<x^2> = 0.496   (kT = 0.5)
    gamma =  10.0: m<v^2> = 0.500, k<x^2> = 0.497   (kT = 0.5)




.. GENERATED FROM PYTHON SOURCE LINES 57-62

Breaking the balance
--------------------
Keep the friction at :math:`\gamma = 1` but give the kicks only half the
fluctuation-dissipation variance, :math:`\gamma k_BT` instead of
:math:`2\gamma k_BT`. The particles cool to half the bath temperature.

.. GENERATED FROM PYTHON SOURCE LINES 62-125

.. code-block:: Python



    @njit
    def trap_drift(s, t, p):
        n = s.shape[0] // 2
        out = np.empty_like(s)
        out[:n] = s[n:]
        out[n:] = (-p[0] * s[:n] - p[1] * s[n:]) / p[2]
        return out


    @njit
    def weak_noise(s, t, p):
        n = s.shape[0] // 2
        out = np.zeros_like(s)
        out[n:] = np.sqrt(p[3]) / p[2]
        return out


    n, dt, chunk = 4000, 0.002, 100
    gamma = 1.0
    weak_params = np.array([k, gamma, m, gamma * kT])
    rng = np.random.default_rng(1908)
    state = np.concatenate([np.full(n, 1.0), np.zeros(n)])
    t_w, s_w = [0.0], [state]
    for i in range(400):  # 400 chunks of 100 steps: t = 0 .. 80
        dW = wiener_increments(chunk, 2 * n, dt, rng=rng)
        _, traj = euler_maruyama_integrate(trap_drift, weak_noise, state, i * chunk * dt, dt, dW, weak_params, chunk)
        state = traj[-1]
        t_w.append((i + 1) * chunk * dt)
        s_w.append(state)
    t_w, s_w = np.array(t_w), np.array(s_w)
    T_weak = m * np.mean(s_w[t_w > 60, n:] ** 2)
    print(f"half-strength noise: m<v^2> = {T_weak:.3f}, not kT = {kT}")

    fig, axes = plt.subplots(1, 3, figsize=(15, 4.3))
    colors = {0.1: "tab:blue", 1.0: "tab:orange", 10.0: "tab:green"}
    for gamma, (t, x, v) in results.items():
        axes[0].plot(t, x[:, 0], lw=0.7, color=colors[gamma], label=rf"$\gamma = {gamma}$")
        axes[1].plot(t, m * np.mean(v**2, axis=1), color=colors[gamma], label=rf"$\gamma = {gamma}$")
    axes[1].plot(t_w, m * np.mean(s_w[:, n:] ** 2, axis=1), color="gray", label="half-strength noise")
    axes[1].axhline(kT, color="k", ls="--", lw=1)
    axes[1].axhline(kT / 2, color="gray", ls=":", lw=1)
    axes[0].set_xlabel("t")
    axes[0].set_ylabel("x")
    axes[0].set_title("one particle per friction")
    axes[0].legend(fontsize=8)
    axes[1].set_xlabel("t")
    axes[1].set_ylabel(r"$m\langle v^2\rangle$")
    axes[1].set_title("kinetic temperature")
    axes[1].legend(fontsize=8)

    t, x, v = results[1.0]
    late = t > 60
    xs = np.linspace(-1.5, 1.5, 200)
    axes[2].hist(x[late].ravel(), bins=80, density=True, alpha=0.6, label="positions, t > 60")
    axes[2].plot(xs, np.sqrt(k / (2 * np.pi * kT)) * np.exp(-k * xs**2 / (2 * kT)), "k--", label=r"Boltzmann $e^{-kx^2/2k_BT}$")
    axes[2].set_xlabel("x")
    axes[2].set_title(r"equilibrium in the trap ($\gamma = 1$)")
    axes[2].legend(fontsize=8)
    fig.tight_layout()
    plt.show()




.. image-sg:: /api/gallery/statphys/stochastic_processes/images/sphx_glr_plot_langevin_fluctuation_dissipation_001.png
   :alt: one particle per friction, kinetic temperature, equilibrium in the trap ($\gamma = 1$)
   :srcset: /api/gallery/statphys/stochastic_processes/images/sphx_glr_plot_langevin_fluctuation_dissipation_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none

    half-strength noise: m<v^2> = 0.253, not kT = 0.5




.. GENERATED FROM PYTHON SOURCE LINES 126-130

Check
-----
With noise strength 2 gamma kT every friction thermalizes to equipartition,
m<v^2> = k<x^2> = kT; halving the noise halves the temperature.

.. GENERATED FROM PYTHON SOURCE LINES 130-133

.. code-block:: Python

    for kinetic, potential in equipartition.values():
        assert abs(kinetic / kT - 1) < 0.03 and abs(potential / kT - 1) < 0.03
    assert abs(T_weak / (kT / 2) - 1) < 0.03








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

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


.. _sphx_glr_download_api_gallery_statphys_stochastic_processes_plot_langevin_fluctuation_dissipation.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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