.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/van_der_pol_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_van_der_pol_demo.py: Van der Pol oscillator — underdamped inference =============================================== Infer the force field of a 1D Van der Pol oscillator from position-only data. Velocities are reconstructed internally by the ULI estimators. This example demonstrates underdamped Langevin inference: the system has inertia, but only positions are observed — as is typical for experimental tracking data. .. rubric:: Tags synthetic · underdamped · nonlinear · 1D .. GENERATED FROM PYTHON SOURCE LINES 16-18 .. code-block:: Python :dedent: 1 .. GENERATED FROM PYTHON SOURCE LINES 39-51 System definition ------------------ The Van der Pol oscillator: .. math:: \ddot{x} = \mu (1 - x^2)\,\dot{x} - x + \sqrt{2D}\,\xi(t) with :math:`\mu = 3` (strongly nonlinear) and :math:`D_0 = 0.05`. We use a time step :math:`\Delta t = 5 \times 10^{-3}` and a long trajectory to capture ~10 complete oscillation cycles. .. GENERATED FROM PYTHON SOURCE LINES 51-83 .. code-block:: Python from SFI.bases import named_scalar, unit_axes, v_components, x_components from SFI.langevin import UnderdampedProcess from SFI.utils import plotting mu = 3.0 D0 = 0.05 dt = 0.02 Nsteps = 2000 oversampling = 8 # fine substeps → accurate Verlet integration of the SDE seed_data = 0 # Ground-truth force written compositionally: # F(x, v) = mu * (1 - x²) * v - x # Named scalar carries the reference value as a default, so the # simulator binds the parameter automatically. (_x,) = x_components(1) (_v,) = v_components(1) (_ex,) = unit_axes(1) mu_p = named_scalar("mu", default=mu) F_exact = (mu_p * (1 - _x * _x) * _v - _x) * _ex key = random.PRNGKey(seed_data) proc = UnderdampedProcess(F_exact, D=jnp.array([[D0]])) proc.initialize(jnp.array([1.0]), v0=jnp.array([0.0])) key, sub = random.split(key) coll = proc.simulate(dt=dt, Nsteps=Nsteps, key=sub, prerun=200, oversampling=oversampling) t, X, _ = coll.to_arrays(dataset=0) # t:(T,), X:(T, N, d) x = X[:, 0, 0] .. GENERATED FROM PYTHON SOURCE LINES 84-90 Phase portrait --------------- The Van der Pol limit cycle in the :math:`(x, \dot{x})` plane. Velocities are estimated via finite differences for visualisation; the inference engine uses more sophisticated ULI reconstructions. .. GENERATED FROM PYTHON SOURCE LINES 90-91 .. code-block:: Python :dedent: 1 .. image-sg:: /gallery/images/sphx_glr_van_der_pol_demo_001.png :alt: Van der Pol oscillator — trajectory data, Phase portrait (time-coloured), Position time series :srcset: /gallery/images/sphx_glr_van_der_pol_demo_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 111-117 Force inference ---------------- Monomials in :math:`(x, v)` up to degree 3 capture the cubic nonlinearity. :class:`~SFI.inference.UnderdampedLangevinInference` handles velocity reconstruction, diffusion estimation, and force projection. .. GENERATED FROM PYTHON SOURCE LINES 117-134 .. code-block:: Python from SFI.bases import monomials_up_to from SFI import UnderdampedLangevinInference degree = 3 B = monomials_up_to(order=degree, dim=1, include_v=True, rank='vector') inf = UnderdampedLangevinInference(coll) inf.compute_diffusion_constant() # Default preset="auto": this trajectory is clean (Lambda_trace < 0), so the # auto switch selects the sharper "clean"/WeakNoise estimator on its own. inf.infer_force_linear(B) inf.compute_force_error() inf.compare_to_exact(model_exact=proc, data_exact=coll, maxpoints=2000) inf.print_report() .. rst-class:: sphx-glr-script-out .. code-block:: none --- StochasticForceInference Report --- Average diffusion tensor: [[0.05549547]] Measurement noise tensor: [[-4.8412693e-07]] Force estimated information: 4455.1044921875 Force: estimated normalized mean squared error (sampling only): 0.001122307275636466 Normalized MSE (force): 0.0229 Normalized MSE (diffusion): 0.0098 Force Coefficient Table ───────────────────────────────────────────────────────────────────── # Label Coefficient Std.Err SNR Sig ───────────────────────────────────────────────────────────────────── 0 1·e0 -2.44047e-01 1.78127e-01 1.4 · 1 x0·e0 -9.87828e-01 1.56945e-01 6.3 * 2 v0·e0 4.52392e+00 1.58589e-01 28.5 ** 3 x0^2·e0 1.18769e-01 5.85006e-02 2.0 * 4 (x0·v0)·e0 -4.98457e-02 6.97995e-02 0.7 · 5 v0^2·e0 1.80063e-02 2.03494e-02 0.9 · 6 x0^3·e0 -1.16508e-01 4.68651e-02 2.5 * 7 (x0^2·v0)·e0 -4.26904e+00 1.00399e-01 42.5 ** 8 (x0·v0^2)·e0 5.84279e-01 5.35458e-02 10.9 ** 9 v0^3·e0 -1.27913e-01 1.23059e-02 10.4 ** ───────────────────────────────────────────────────────────────────── 10/10 basis functions in support, sig: 7* / 4** / 0*** (|SNR| ≥ 2 / 10 / 100) (Std.err. reflects sampling error only; discretization bias is not included.) .. GENERATED FROM PYTHON SOURCE LINES 135-140 Force validation ----------------- Scatter plot of exact vs inferred force, evaluated along the ULI-reconstructed trajectory. .. GENERATED FROM PYTHON SOURCE LINES 140-143 .. code-block:: Python F_true, F_pred = inf.force_comparison_arrays(model_exact=proc) .. image-sg:: /gallery/images/sphx_glr_van_der_pol_demo_002.png :alt: Underdamped force inference — validation, $r=0.993$ MSE=0.012, Force along trajectory :srcset: /gallery/images/sphx_glr_van_der_pol_demo_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 168-173 Bootstrapped trajectory ------------------------- Simulating from the inferred model and comparing the phase portrait provides a generative validation. .. GENERATED FROM PYTHON SOURCE LINES 173-182 .. code-block:: Python key_boot = random.PRNGKey(1) coll_boot, _ = inf.simulate_bootstrapped_trajectory( key_boot, oversampling=oversampling, simulate=True ) _, X_boot, _ = coll_boot.to_arrays(dataset=0) x_boot = X_boot[:, 0, 0] v_fd_boot = coll_boot.velocity_array(dataset=0, scheme="central")[:, 0, 0] .. image-sg:: /gallery/images/sphx_glr_van_der_pol_demo_003.png :alt: Generative validation — bootstrapped trajectory, Phase portrait: data vs bootstrap, Time series comparison :srcset: /gallery/images/sphx_glr_van_der_pol_demo_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 208-219 Diagnostics ----------- The **underdamped** residual compares the symmetric finite-difference acceleration :math:`\hat a = (x_{t+1} - 2x_t + x_{t-1})/\Delta t^2` (the same kinematic the force estimator fits) to the inferred force :math:`\hat F(\hat x, \hat v)`, standardised by the noise scale of a second difference of integrated white noise, :math:`\tfrac23\,(2\hat D)/\Delta t`. Adjacent accelerations overlap, so every second sample is used to keep the residuals independent. A well-specified cubic-in-:math:`(x,v)` model then passes all checks. .. GENERATED FROM PYTHON SOURCE LINES 219-225 .. code-block:: Python from SFI.diagnostics import assess, plot_summary report = assess(inf, level="standard") report.print_summary() .. image-sg:: /gallery/images/sphx_glr_van_der_pol_demo_004.png :alt: Van der Pol — diagnostics (underdamped linear model), Q--Q plot, Residual histogram, Residual ACF :srcset: /gallery/images/sphx_glr_van_der_pol_demo_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none === SFI diagnostics report === backend : linear regime : UD n_obs : 999 n_particles: 1 d: 1 level : standard -- Residuals -- mean = +0.0221 std = 1.0195 skew = -0.071 kurt-3 = -0.129 (n=999) ✓ normality ks stat=0.0294 p=0.347 ✗ autocorr ljung_box stat=62.9 p=2.54e-06 ✗ autocorr ljung_box_squared stat=46.1 p=0.00078 predicted NMSE = 0.00112 realised NMSE = 0.00534 χ² z = +0.87 (|z|>5 ⇒ bias) -- Flags -- ! [autocorr/ljung_box] p=2.54e-06 < 0.01 — missing time-correlated feature — widen the basis; if it persists, suspect coarse sampling: the parametric estimator (infer_force) extends the usable Δt ! [autocorr/ljung_box_squared] p=7.80e-04 < 0.01 — diffusion mis-estimated or state-dependent — try the other compute_diffusion_constant method or a state-dependent diffusion basis; the parametric estimators profile (D, Λ) .. GENERATED FROM PYTHON SOURCE LINES 231-233 Thumbnail --------- .. GENERATED FROM PYTHON SOURCE LINES 233-236 .. code-block:: Python stamp_output() .. image-sg:: /gallery/images/sphx_glr_van_der_pol_demo_005.png :alt: van der pol demo :srcset: /gallery/images/sphx_glr_van_der_pol_demo_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [Generated: 2026-06-30 10:05] .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 17.277 seconds) .. rst-class:: sphx-glr-example-tags 🏷 Tags: synthetic, underdamped, nonlinear, 1D .. _sphx_glr_download_gallery_van_der_pol_demo.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: van_der_pol_demo.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: van_der_pol_demo.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: van_der_pol_demo.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_