.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/entropy_production_underdamped_demo.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_gallery_entropy_production_underdamped_demo.py: Inertial dissipation: entropy production from positions only ============================================================= The :doc:`overdamped entropy-production demo ` measured the dissipation of a driven optical trap. Here the same bead keeps its **inertia**: the state is :math:`(\mathbf{x}, \mathbf{v})`, the bath exerts friction and noise on the velocity, .. math:: \frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t} = \mathbf{v}, \qquad \frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t} = \underbrace{-k\,\mathbf{x} + \varepsilon\,(\hat z \times \mathbf{x})}_{\text{trap + curl drive}} \;\underbrace{-\;\gamma\,\mathbf{v}}_{\text{friction}} \;+\; \sqrt{2D}\,\boldsymbol{\xi}(t), and — crucially — **only positions are observed**. Underdamped Langevin Inference reconstructs the velocities, fits the full acceleration field :math:`F(\mathbf{x},\mathbf{v})`, and splits it by time-reversal parity (:math:`\mathbf{x}` even, :math:`\mathbf{v}` odd): a *reversible* part :math:`F^{+}` (trap and drive) and an *irreversible* part :math:`F^{-}` (friction). The entropy production is then the Stratonovich log path-probability ratio of the fitted model, evaluated on the data. Inertia is not a detail here. For this linear system the exact rate is .. math:: \sigma(\gamma) \;=\; \frac{2\,\varepsilon^2\,\gamma}{\gamma^2 k - \varepsilon^2}, which *grows* as friction decreases — diverging at the stability boundary :math:`\gamma \to \varepsilon/\sqrt{k}` where the curl overcomes the trap — and only settles onto the overdamped result :math:`2\varepsilon^2/(\gamma k)` at strong friction. Analysing inertial data with an overdamped model misses this amplification. Across this page we will: #. simulate the inertial driven trap, recording positions only; #. fit :math:`F(\mathbf{x},\mathbf{v})` with ULI and let PASTIS find the trap, curl and friction terms; #. split the fit by time-reversal parity with :meth:`~SFI.inference.UnderdampedLangevinInference.time_reversal_split`; #. **measure the entropy production** with :meth:`~SFI.inference.UnderdampedLangevinInference.compute_entropy_production`, cross-check it against the simulator's ground truth, and remove the same-sample plug-in bias by cross-fitting; #. sweep the friction: inertia amplifies dissipation, and the estimator tracks the exact :math:`\sigma(\gamma)`; #. push the sampling interval past the velocity correlation time and watch the reconstruction fail — the :math:`\Delta t \lesssim \tau_v` prerequisite. .. rubric:: Tags synthetic · underdamped · linear · 2D · non-equilibrium · thermodynamics .. GENERATED FROM PYTHON SOURCE LINES 62-65 .. code-block:: Python :dedent: 1 .. GENERATED FROM PYTHON SOURCE LINES 87-102 An inertial bead in the driven trap ----------------------------------- Same trap, same faint rotational drive as the overdamped page — but the bead now carries momentum (think a levitated microsphere in low-pressure gas rather than a colloid in water). We pick a moderately underdamped regime, :math:`\gamma = \sqrt{k}` (quality factor near one), where the exact rate :math:`\sigma = 2\varepsilon^2\gamma/(\gamma^2 k - \varepsilon^2)` is already **1.33×** the overdamped formula. The simulator integrates the phase-space dynamics but the returned collection contains **positions only** — exactly what a camera sees. ``compute_observables=True`` makes it accumulate the ground-truth information and entropy along the way, using the true internal velocities (which are otherwise discarded). .. GENERATED FROM PYTHON SOURCE LINES 102-148 .. code-block:: Python import SFI from SFI.bases import unit_axes, v_components, x_components from SFI.langevin import UnderdampedProcess k_true = 1.0 # trap stiffness eps_true = 0.5 # rotational drive: F_curl = eps * (z_hat x x) gamma_true = 1.0 # friction D_true = 0.5 # velocity-space diffusion dt = 0.01 # time between recorded frames Nsteps = 320_000 # frames (3200 velocity relaxation times) x0c, x1c = x_components(2) v0c, v1c = v_components(2) e0, e1 = unit_axes(2) def make_trap(eps, gamma): """Trap + curl drive + friction, as a basis expression F(x, v).""" return ( (-k_true * x0c - eps * x1c - gamma * v0c) * e0 + (-k_true * x1c + eps * x0c - gamma * v1c) * e1 ) def simulate_trap(eps, gamma, seed, *, dt=dt, Nsteps=Nsteps, oversampling=8): proc = UnderdampedProcess(make_trap(eps, gamma), D=D_true) proc.initialize(jnp.array([0.5, 0.0], dtype=jnp.float32)) coll_run = proc.simulate( dt=dt, Nsteps=Nsteps, key=random.PRNGKey(seed), prerun=500, oversampling=oversampling, compute_observables=True, ) return coll_run, proc def sigma_exact(eps, gamma): return 2.0 * eps**2 * gamma / (gamma**2 * k_true - eps**2) coll, proc_true = simulate_trap(eps_true, gamma_true, seed=3) print(f"Simulated {coll.T} frames at dt = {dt} (positions only)") print(f"Exact entropy production rate: sigma = {sigma_exact(eps_true, gamma_true):.4f} kB / time") print(f"Overdamped formula would give: {2 * eps_true**2 / (gamma_true * k_true):.4f} kB / time") .. image-sg:: /gallery/images/sphx_glr_entropy_production_underdamped_demo_001.png :alt: Position vs time: smoother than overdamped, The bead's path: inertial loops :srcset: /gallery/images/sphx_glr_entropy_production_underdamped_demo_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Simulated 320000 frames at dt = 0.01 (positions only) Exact entropy production rate: sigma = 0.6667 kB / time Overdamped formula would give: 0.5000 kB / time .. GENERATED FROM PYTHON SOURCE LINES 170-178 Fit the acceleration field from positions ----------------------------------------- The workflow mirrors the overdamped one — the engine transparently handles the velocity reconstruction. We hand ULI an over-complete polynomial basis in :math:`(\mathbf{x}, \mathbf{v})` and let PASTIS decide: it should keep the linear trap terms, the faint curl, and the friction — and nothing else. .. GENERATED FROM PYTHON SOURCE LINES 178-196 .. code-block:: Python from SFI.bases import monomials_up_to inf = SFI.UnderdampedLangevinInference(coll) inf.compute_diffusion_constant() B = monomials_up_to(2, dim=2, include_constant=True, include_v=True, rank="vector") inf.infer_force_linear(B) inf.compute_force_error() coeffs_full = np.asarray(inf.force_coefficients_full) stderr_full = np.asarray(inf.force_coefficients_stderr) inf.sparsify_force(criterion="PASTIS") inf.compute_force_error() inf.print_report() .. image-sg:: /gallery/images/sphx_glr_entropy_production_underdamped_demo_002.png :alt: 30 candidate terms: trap, curl and friction recovered :srcset: /gallery/images/sphx_glr_entropy_production_underdamped_demo_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none --- StochasticForceInference Report --- Average diffusion tensor: [[ 5.0369191e-01 -3.4161386e-04] [-3.4161386e-04 5.0351131e-01]] Measurement noise tensor: [[ 5.1027321e-10 -6.0447689e-12] [-6.0447689e-12 3.9476983e-10]] Force estimated information: 3711.53662109375 Force: estimated normalized mean squared error (sampling only): 0.000808290684722677 Force Coefficient Table ─────────────────────────────────────────────────────────────────── # Label Coefficient Std.Err SNR Sig ─────────────────────────────────────────────────────────────────── 2 x0·e0 -1.00016e+00 2.16500e-02 46.2 ** 3 x0·e1 5.18507e-01 2.49182e-02 20.8 ** 4 x1·e0 -4.84258e-01 2.51104e-02 19.3 ** 5 x1·e1 -9.95782e-01 2.18199e-02 45.6 ** 6 v0·e0 -9.95473e-01 2.50260e-02 39.8 ** 9 v1·e1 -1.03102e+00 2.52615e-02 40.8 ** ─────────────────────────────────────────────────────────────────── 6/30 basis functions in support, sig: 6* / 6** / 0*** (|SNR| ≥ 2 / 10 / 100) (Std.err. reflects sampling error only; discretization bias is not included.) Zeroed (24): 1·e0, 1·e1, v0·e1, v1·e0, x0^2·e0, x0^2·e1, (x0·x1)·e0, (x0·x1)·e1, (x0·v0)·e0, (x0·v0)·e1, (x0·v1)·e0, (x0·v1)·e1, x1^2·e0, x1^2·e1, (x1·v0)·e0, (x1·v0)·e1, (x1·v1)·e0, (x1·v1)·e1, v0^2·e0, v0^2·e1, (v0·v1)·e0, (v0·v1)·e1, v1^2·e0, v1^2·e1 .. GENERATED FROM PYTHON SOURCE LINES 217-226 Split the dynamics by the arrow of time --------------------------------------- Under time reversal positions are even and velocities odd, so the fitted field splits into :math:`\hat F^{\pm}(\mathbf{x},\mathbf{v}) = \tfrac12[\hat F(\mathbf{x},\mathbf{v}) \pm \hat F(\mathbf{x},-\mathbf{v})]`. For this system :math:`\hat F^{+}` is the trap plus the curl drive (all positional) and :math:`\hat F^{-}` is the friction — the part that knows which way time flows. .. GENERATED FROM PYTHON SOURCE LINES 226-236 .. code-block:: Python F_even, F_odd = inf.time_reversal_split() x_probe = jnp.array([[1.0, 0.0]]) v_probe = jnp.array([[0.0, 1.0]]) print("F (x=(1,0), v=(0,1)) =", np.asarray(inf.force_inferred(x_probe, v=v_probe))[0]) print("F_even (reversible part) =", np.asarray(F_even(x_probe, v_probe))[0]) print("F_odd (irreversible part)=", np.asarray(F_odd(x_probe, v_probe))[0]) .. image-sg:: /gallery/images/sphx_glr_entropy_production_underdamped_demo_003.png :alt: $\hat F^{+}(\mathbf{x}, 0)$: trap + drive (even), $\hat F^{-}$: friction (odd) :srcset: /gallery/images/sphx_glr_entropy_production_underdamped_demo_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none F (x=(1,0), v=(0,1)) = [-1.0001588 -0.5125175] F_even (reversible part) = [-1.0001588 0.51850665] F_odd (irreversible part)= [ 0. -1.0310241] .. GENERATED FROM PYTHON SOURCE LINES 264-280 Measure the dissipation — from positions only --------------------------------------------- :meth:`~SFI.inference.UnderdampedLangevinInference.compute_entropy_production` evaluates :math:`\Delta\hat S = \sum_t (\mathrm{d}\hat v_t \circ - \hat F^{+}\mathrm{d}t)^\top \bar D^{-1} \hat F^{-}` with the reconstructed kinematics. Two subtleties, both reported by the method: - friction makes the odd sector *informative but barely irreversible*: the log ratio's contributions largely cancel, so the error bar is set by the odd-sector fluctuation scale (``Q_odd``), estimated analytically and by block variance; - evaluating the functional with coefficients fitted **on the same trajectory** adds a small positive :math:`O(1/\tau_N)` plug-in bias. **Cross-fitting** — fit on one half, evaluate on the other — removes it, via the ``coefficients`` argument. .. GENERATED FROM PYTHON SOURCE LINES 280-317 .. code-block:: Python out = inf.compute_entropy_production() S_sim = coll.datasets[0].meta["observables"]["entropy"] tauN = out["tauN"] sigma = sigma_exact(eps_true, gamma_true) # Cross-fitting: two halves, two fits, evaluate each on the other half. from SFI import TrajectoryCollection X_all = jnp.asarray(X_ud) T_half = X_all.shape[0] // 2 halves = [ TrajectoryCollection.from_arrays(X=X_all[:T_half], dt=dt), TrajectoryCollection.from_arrays(X=X_all[T_half:], dt=dt), ] fits = [] for coll_h in halves: inf_h = SFI.UnderdampedLangevinInference(coll_h) inf_h.compute_diffusion_constant() inf_h.infer_force_linear(B) fits.append(inf_h) out_xf = [ fits[1].compute_entropy_production(coefficients=fits[0].force_coefficients_full), fits[0].compute_entropy_production(coefficients=fits[1].force_coefficients_full), ] Sdot_crossfit = 0.5 * (out_xf[0]["Sdot"] + out_xf[1]["Sdot"]) err_crossfit = 0.5 * (out_xf[0]["Sdot_error"] ** 2 + out_xf[1]["Sdot_error"] ** 2) ** 0.5 print(f"exact rate : {sigma:.4f} kB / time") print(f"inferred (same-sample) : {out['Sdot']:.4f} +- {out['Sdot_error']:.4f}" f" [odd dof bias {out['Sdot_bias']:.4f}]") print(f"inferred (cross-fitted) : {Sdot_crossfit:.4f} +- {err_crossfit:.4f}") print(f"true model on this trajectory : {S_sim / tauN:.4f}") print(f"irreversible fraction of the odd sector: {out['entropy_odd_ratio']:.3f}") .. image-sg:: /gallery/images/sphx_glr_entropy_production_underdamped_demo_004.png :alt: Dissipation from positions only :srcset: /gallery/images/sphx_glr_entropy_production_underdamped_demo_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none exact rate : 0.6667 kB / time inferred (same-sample) : 0.6624 +- 0.0423 [odd dof bias 0.0013] inferred (cross-fitted) : 0.6785 +- 0.0517 true model on this trajectory : 0.6562 irreversible fraction of the odd sector: 0.246 .. GENERATED FROM PYTHON SOURCE LINES 336-346 Inertia amplifies dissipation ----------------------------- Now sweep the friction at fixed drive. Watch the exact rate: far from the overdamped intuition that less friction means less dissipation, :math:`\sigma(\gamma)` **grows** as :math:`\gamma` decreases — the freely-circulating bead extracts more work from the drive — and diverges at the stability boundary :math:`\gamma = \varepsilon/\sqrt{k}`. The estimator, still fed positions only, tracks the exact curve; the overdamped formula (dashed) is off by 30% already at :math:`\gamma = 1`. .. GENERATED FROM PYTHON SOURCE LINES 346-364 .. code-block:: Python B_lin = monomials_up_to(1, dim=2, include_constant=False, include_v=True, rank="vector") def entropy_of(coll_run): inf_run = SFI.UnderdampedLangevinInference(coll_run) inf_run.compute_diffusion_constant() inf_run.infer_force_linear(B_lin) return inf_run.compute_entropy_production() gamma_sweep = [0.75, 1.0, 1.5, 2.0, 3.0, 5.0] sweep = {} for i, gam in enumerate(gamma_sweep): coll_g, _ = simulate_trap(eps_true, gam, seed=20 + i, Nsteps=80_000) sweep[gam] = entropy_of(coll_g) .. image-sg:: /gallery/images/sphx_glr_entropy_production_underdamped_demo_005.png :alt: Inertia amplifies dissipation :srcset: /gallery/images/sphx_glr_entropy_production_underdamped_demo_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 389-399 When sampling is too coarse --------------------------- Everything above relied on resolving the velocity correlation time :math:`\tau_v = 1/\gamma`. Subsampling the same recording stretches the effective :math:`\Delta t`; beyond :math:`\gamma\Delta t \approx 0.5` the secant velocities stop resembling velocities, the estimator (which warns at exactly that threshold) degrades, and the measurement slides toward a coarse-grained lower bound. There is no free lunch: **inertial dissipation is only measurable when inertia is resolved.** .. GENERATED FROM PYTHON SOURCE LINES 399-410 .. code-block:: Python strides = [1, 5, 10, 20, 40, 80] resolution = [] for s in strides: coll_s = TrajectoryCollection.from_arrays(X=X_all[::s], dt=dt * s) resolution.append((gamma_true * dt * s, entropy_of(coll_s))) for gdt, o in resolution: print(f"gamma*dt = {gdt:5.2f}: Sdot = {o['Sdot']:7.4f} +- {o['Sdot_error']:.4f}") .. image-sg:: /gallery/images/sphx_glr_entropy_production_underdamped_demo_006.png :alt: Velocities must be resolved: $\Delta t \lesssim \tau_v$ :srcset: /gallery/images/sphx_glr_entropy_production_underdamped_demo_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none gamma*dt = 0.01: Sdot = 0.6632 +- 0.0425 gamma*dt = 0.05: Sdot = 0.7219 +- 0.0466 gamma*dt = 0.10: Sdot = 0.7975 +- 0.0517 gamma*dt = 0.20: Sdot = 0.9738 +- 0.0630 gamma*dt = 0.40: Sdot = 1.3950 +- 0.0893 gamma*dt = 0.80: Sdot = 2.0682 +- 0.1287 .. GENERATED FROM PYTHON SOURCE LINES 429-443 Next steps ---------- - **The overdamped companion page**: the same trap without inertia, the near-equilibrium detection limit, and the parametric energy fit — :doc:`/gallery/entropy_production_demo`. - **Theory**: the odd-parity split, the estimator's provisional error bars and the cross-fitting recipe — :doc:`/inference/stochastic_thermodynamics`. - **Driven or inertial?** The dynamics-order classifier tells a rotational current from momentum — :doc:`/gallery/dynamics_order_demo`. - **Underdamped inference in depth**: presets, kinematics and diagnostics — :doc:`/inference/underdamped`. .. GENERATED FROM PYTHON SOURCE LINES 445-447 Thumbnail --------- .. GENERATED FROM PYTHON SOURCE LINES 447-450 .. code-block:: Python stamp_output() .. image-sg:: /gallery/images/sphx_glr_entropy_production_underdamped_demo_007.png :alt: entropy production underdamped demo :srcset: /gallery/images/sphx_glr_entropy_production_underdamped_demo_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [Generated: 2026-07-03 11:13] .. rst-class:: sphx-glr-timing **Total running time of the script:** (3 minutes 57.575 seconds) .. rst-class:: sphx-glr-example-tags 🏷 Tags: synthetic, underdamped, linear, 2D, non-equilibrium, thermodynamics .. _sphx_glr_download_gallery_entropy_production_underdamped_demo.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: entropy_production_underdamped_demo.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: entropy_production_underdamped_demo.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: entropy_production_underdamped_demo.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_