TinyUDE: Solver-Free Universal Differential Equations on Microcontrollers via Lie-Taylor Jet Matching

# TinyUDE: Solver-Free Universal Differential Equations on Microcontrollers via Lie–Taylor Jet Matching
Source: [https://arxiv.org/html/2609.26972](https://arxiv.org/html/2609.26972)
Pranavanath BalamuraliAffiliation:Department of Electrical and Computer EngineeringAffiliation:The University of Texas at AustinEmail:[pb27677@utexas\.edu](mailto:)Hrishi KamireddyAffiliation:Department of Computer ScienceAffiliation:The University of Texas at AustinEmail:[hk27965@utexas\.edu](mailto:)

Abstract

Training Universal Differential Equations \(UDEs\) traditionally relies on backpropagating through numerical ODE solvers, creating memory footprints far exceeding the capabilities of edge microcontrollers\. We present Lie–Taylor jet matching, a solver\-free training framework that fits a hybrid vector field directly to the first and second time\-derivatives of observed system states\. These derivatives, the truncated Lie–Taylor jet, are estimated online via Savitzky–Golay filtering, yielding fully analytic gradients without automatic differentiation software\. We evaluate whether eliminating the solver compromises accuracy against a conventional baseline \(fixed\-step RK4 integration, multiple shooting, exact discrete adjoints, Adam\) sharing identical dynamics, noise models, network architectures, and metrics\. While naive derivative matching degrades under sensor noise, our noise\-adaptive mechanisms close and reverse this gap: full\-rate phase\-shifted sampling, a reservoir buffer, cosine\-annealed optimization with weight averaging, on\-device noise estimation, and polynomial\-misfit quality gating\. On a damped pendulum and chaotic double pendulum, our method matches or exceeds baseline accuracy at matched data windows and recovers unmodeled damping coefficients\. Across noise levels from0%0\\%to5%5\\%, it attains a geometric\-mean relative field error of0\.65×0\.65\\timesthat of the baseline within108108kB of static memory, compared with megabytes of solver tape\. On an ESP32 microcontroller, the on\-device run reaches a field error of0\.00200\.0020and recovers the damping coefficient toc^=0\.400\\hat\{c\}=0\.400\(true0\.4000\.400\) within61\.361\.3kB of static memory and7\.247\.24ms per update \(18\.1%18\.1\\%duty cycle at2525Hz\), confirming real\-time on\-device training is feasible without a numerical solver\.

Keywords:Universal Differential Equations⋅\\cdotSolver\-Free Training⋅\\cdotJet Matching⋅\\cdotSavitzky–Golay Filtering⋅\\cdotEmbedded AI⋅\\cdotLie\-Taylor Jet

## 1Introduction

### 1\.1Motivation & Problem Statement

Physical systems operating in field conditions \(robotic manipulators, unmanned aerial vehicles, electric motors, remote sensing nodes\) frequently deviate from initial nominal models\. Mechanical components wear over time, aerodynamic properties change, loads vary, and temperature shifts alter internal friction\. Universal Differential Equations \(UDEs\), differential equations whose right\-hand side is partly a neural network trained inside the equation, provide a powerful framework to address these variations by combining known physics, expressed as an analytic vector fieldfknownf\_\{\\mathrm\{known\}\}, with a neural network residualNθN\_\{\\theta\}designed to capture unmodeled dynamics\[[1](https://arxiv.org/html/2609.26972#bib.bib1)\]\. This hybrid formulation preserves physical interpretability and structural extrapolation while granting the model sufficient capacity to capture unexpected disturbances, drag, and dissipation\.

Despite their advantages, the standard training pipeline for UDEs creates significant computational hurdles for edge deployments\. Typically, the hybrid system is integrated inside an optimization loop, and gradients are computed by backpropagating through the differential equation solver using unrolled automatic differentiation or continuous adjoint states \(the adjoint method obtains parameter gradients by integrating a second, backward\-in\-time differential equation for the loss sensitivities\)\[[2](https://arxiv.org/html/2609.26972#bib.bib2)\]\. This is often stabilized over long time horizons using multiple shooting\[[3](https://arxiv.org/html/2609.26972#bib.bib5)\], which splits a long trajectory into short windows that are each re\-initialized from data so the solver never integrates far from the observations\. This standard workflow assumes workstation\-class computational resources, requiring dedicated automatic differentiation libraries, megabytes of intermediate state storage \(tapes\), and full offline trajectory access\.

In this work, we investigate whether Ordinary Differential Equation \(ODE\) solvers and intermediate backpropagation tapes can be eliminated entirely\. By matching estimated time\-derivatives in a streaming framework, we examine how solver\-free learning compares to traditional adjoint methods in terms of model accuracy\. We demonstrate that when derivative estimates are processed with appropriate statistical filtering, direct derivative matching eliminates the collocation accuracy penalty\. Specifically, noise\-adaptive Lie–Taylor jet matching outperforms validated adjoint baselines while reducing training memory requirements by roughly two orders of magnitude\. Because direct numerical differentiation inherently amplifies measurement noise as𝒪⁡\(σ/Δ​t\)\\mathcal\{O\}\(\\sigma/\\Delta t\), standard least\-squares regression on noisy targets produces biased estimates\. The mechanisms introduced in Section[3](https://arxiv.org/html/2609.26972#S3)address this by creating statistically independent derivative targets, averaging over streaming samples, operating filters at noise\-optimal parameters, and rejecting corrupted data points\.

### 1\.2Edge Constraints: Memory Limitations and Solver Fragility

Two primary obstacles restrict traditional UDE training on resource\-constrained microcontrollers:

1. 1\.Memory Footprint of Backpropagation\.Differentiating through aKK\-step numerical solver requires storing intermediate solver stages and network activations for every step in the forward pass\. In our baseline experiments \(Section[4](https://arxiv.org/html/2609.26972#S4)\), a modest network with 354 parameters trained using a multiple shooting window ofK=40K=40requires5\.95\.9MB of memory tape and peaks at9\.19\.1MB of heap usage\. Extending the window toK=80K=80to handle higher noise increases memory usage to11\.811\.8MB\. Standard microcontrollers, which typically possess around320320kB of static random\-access memory \(SRAM\), cannot host these operations\. Consequently, existing embedded ML frameworks focus almost exclusively on model inference\[[4](https://arxiv.org/html/2609.26972#bib.bib14),[5](https://arxiv.org/html/2609.26972#bib.bib15)\], while on\-device training methods are generally restricted to fine\-tuning final classification layers rather than updating dynamical systems\[[6](https://arxiv.org/html/2609.26972#bib.bib16)\]\.
2. 2\.Numerical Integrator Divergence\.Embedded numerical integration inside an optimization loop tightly couples model parameter updates to solver stability\. Temporary parameter excursions during training can cause intermediate trajectories to diverge, leading to unstable gradient steps\. Furthermore, computing each gradient requires4​K4Kvector\-field evaluations in both forward and reverse passes\. Removing the solver avoids solver instability but introduces a different challenge: training targets must be estimated directly from noisy sensor data\. Without adequate filtering, point\-to\-point collocation fails under noise, and in our experiments, single\-step collocation \(K=1K=1\) degraded to a relative field error of0\.820\.82under5%5\\%noise\. Therefore, solver\-free methods require robust statistical filtering to generate reliable derivative targets\.

### 1\.3Summary of Contributions

We propose noise\-adaptive Lie–Taylor jet matching and present the following key contributions:

- •Solver\-Free Training with Exact Analytic Gradients\.We introduce a solver\-free training approach that fits hybrid vector fields to first\- and second\-order time derivatives\(x˙,x¨\)\(\\dot\{x\},\\ddot\{x\}\)derived from Savitzky–Golay \(SG\) filtering\[[7](https://arxiv.org/html/2609.26972#bib.bib6)\], a sliding\-window least\-squares polynomial fit whose derivatives at the window center serve as smoothed derivative estimates\. The second derivative is modeled using the Lie derivative identityx¨=JFθ​Fθ\\ddot\{x\}=J\_\{F\_\{\\theta\}\}F\_\{\\theta\}, where the Lie derivative is the rate of change of a quantity along the flow of the vector field\[[8](https://arxiv.org/html/2609.26972#bib.bib55)\]\. Both terms yield closed\-form analytic gradients, avoiding the need for automatic differentiation software\. We validate these gradients against complex\-step differentiation \(a technique that recovers derivatives of a real function to machine precision by evaluating it at a small imaginary perturbation, free of the cancellation error of finite differences\) across a 32\-check suite, achieving agreement to machine precision\.
- •On\-Device Noise\-Adaptive Processing\.We implement five lightweight mechanisms designed for microcontrollers to mitigate derivative noise: 1. 1\.*Full\-rate phase sampling*: By sliding the filtering window every few clock ticks while maintaining tap spacingmm, we extract approximately50×50\\timesmore statistically independent derivative targets from raw sensor streams, requiring only a4\.44\.4kB ring buffer\. 2. 2\.*Reservoir buffering*: A20482048–40964096sample buffer \(4848–9696kB\) maintained by reservoir sampling, which keeps a uniform random subset of an unbounded stream in fixed memory, replaces a small rolling minibatch and prevents the network from memorizing the noise on a handful of targets\. 3. 3\.*Averaged Adam optimization*: We use Adam \(adaptive moment estimation, a first\-order optimizer with per\-parameter step sizes\) with a cosine\-annealed learning rate \(a step size that decays along a half cosine so that updates settle\) and Polyak weight averaging implemented as an exponential moving average \(EMA\) of the iterates, which cancels residual optimization noise\[[9](https://arxiv.org/html/2609.26972#bib.bib9),[10](https://arxiv.org/html/2609.26972#bib.bib10)\], all within4\.24\.2kB of state memory\. 4. 4\.*Adaptive filter time\-scaling*: An on\-device noise estimator computes sensor noise directly from fit residuals \(0\.04%0\.04\\%empirical error\) to set optimal filter time scales according toΔ​t∗∝σ1/5/ωc\\Delta t^\{\\ast\}\\propto\\sigma^\{1/5\}/\\omega\_\{c\}, while gating second\-order derivative contributions based on signal\-to\-noise ratios\. 5. 5\.*Polynomial misfit quality gating*: We reject sampling windows where local polynomial assumptions fail, preventing corrupted derivative targets during fast non\-linear state transitions\.
- •Experimental Validation and Performance Evaluation\.We evaluate the approach against an exact discrete adjoint baseline sharing identical plants, noise models, neural networks, and evaluation metrics\. Jet matching achieves a geometric\-mean field error ratio of0\.65×0\.65\\timesrelative to the reference baseline across noise levels from0%0\\%to5%5\\%\(0\.82×0\.82\\timesagainst a baseline configured withK=80K=80\)\. On a chaotic double pendulum, the proposed method yields lower field error \(0\.01160\.0116vs\.0\.01370\.0137\) and accurately recovers unmodeled damping coefficients \(c^=0\.388\\hat\{c\}=0\.388–0\.3990\.399compared to the true value of0\.400\.40\), all within a static memory allocation of108108kB\.

## 2Related Work

### 2\.1Universal Differential Equations and Neural ODEs

Neural ODEs place a neural network on the right\-hand side of a differential equation and train it by differentiating through the numerical solution, either by unrolling the solver steps or by integrating a continuous adjoint equation backward in time\[[2](https://arxiv.org/html/2609.26972#bib.bib2)\]\. The framework has been extended to irregularly sampled time series\[[11](https://arxiv.org/html/2609.26972#bib.bib19)\], augmented with auxiliary dimensions to relax topological constraints on the learned flow\[[12](https://arxiv.org/html/2609.26972#bib.bib20)\], and surveyed comprehensively by a recent study\[[13](https://arxiv.org/html/2609.26972#bib.bib18)\]\. UDEs refine the idea for scientific settings by combining a known mechanistic vector field with a neural term that models only the unknown dynamics\[[1](https://arxiv.org/html/2609.26972#bib.bib1)\], retaining interpretability and extrapolation where the physics is trusted\. A parallel line builds physical structure into the network itself: Hamiltonian neural networks learn a scalar energy whose symplectic gradient gives the dynamics\[[14](https://arxiv.org/html/2609.26972#bib.bib4)\], and Lagrangian formulations do the same through the Euler–Lagrange equations\[[15](https://arxiv.org/html/2609.26972#bib.bib26),[16](https://arxiv.org/html/2609.26972#bib.bib27)\]\. Physics\-informed neural networks \(PINNs\) instead penalize the residual of the governing equations at collocation points using automatic differentiation of a neural representation of the solution\[[17](https://arxiv.org/html/2609.26972#bib.bib3),[18](https://arxiv.org/html/2609.26972#bib.bib28)\]\.

All of these approaches share the training obstacle that motivates this work: the gradient is obtained by reverse\-mode differentiation through a numerical integration, whose memory grows with the number of solver steps\. The literature has attacked the cost from several directions\. Checkpointing and adaptive\-checkpoint adjoints trade recomputation for memory\[[19](https://arxiv.org/html/2609.26972#bib.bib21),[20](https://arxiv.org/html/2609.26972#bib.bib22)\]\. The distinction between discretize\-then\-optimize versus optimize\-then\-discretize clarifies when a continuous adjoint yield exact gradients for the discrete solver used in practice\[[21](https://arxiv.org/html/2609.26972#bib.bib23)\]\. Meanwhile, multiple shooting splits long trajectories into short windows, allowing solver\-based training to tolerate noisy records\[[3](https://arxiv.org/html/2609.26972#bib.bib5),[22](https://arxiv.org/html/2609.26972#bib.bib24)\]\. Collocation\-style training of neural ODEs, which fits the field to smoothed derivative estimates rather than integrating, has been studied as a fast initialization or alternative to adjoint training\[[23](https://arxiv.org/html/2609.26972#bib.bib25)\]\. Our baseline \(Section[4](https://arxiv.org/html/2609.26972#S4)\) is a careful instance of the adjoint recipe: fixed\-step Runge–Kutta integration\[[24](https://arxiv.org/html/2609.26972#bib.bib52),[25](https://arxiv.org/html/2609.26972#bib.bib17)\], multiple shooting, an exact discrete adjoint, and its own shooting\-length sweep quantifies the memory that noise robustness costs\. Jet matching keeps the UDE decomposition but abandons the solve\-inside\-the\-loss construction entirely\. Our chaotic test bed, the damped double pendulum, is a standard benchmark for sensitive dependence on initial conditions\[[26](https://arxiv.org/html/2609.26972#bib.bib54),[27](https://arxiv.org/html/2609.26972#bib.bib53)\]\.

### 2\.2On\-Device Machine Learning and Memory Constraints

Deploying machine learning on microcontroller units \(MCUs\) has matured into a distinct field, TinyML, with interpreter runtimes that execute quantized networks in tens of kilobytes\[[4](https://arxiv.org/html/2609.26972#bib.bib14),[28](https://arxiv.org/html/2609.26972#bib.bib48),[29](https://arxiv.org/html/2609.26972#bib.bib49)\]\. Additionally, hardware\-aware architecture search fits large classifiers into microcontroller flash\[[5](https://arxiv.org/html/2609.26972#bib.bib15)\], and standardized benchmarks for inference latency and energy\[[30](https://arxiv.org/html/2609.26972#bib.bib47)\]\. On\-device training is far less developed\[[31](https://arxiv.org/html/2609.26972#bib.bib51)\]\. Memory\-efficient transfer learning freezes the weights and updates only biases and small adapter modules to avoid storing activations\[[32](https://arxiv.org/html/2609.26972#bib.bib50)\]\. Sparse\-update training under a 256 kB budget selects which layers and channels may change at all\[[6](https://arxiv.org/html/2609.26972#bib.bib16)\], and online learning on microcontrollers has been demonstrated for incremental classifier adaptation\[[28](https://arxiv.org/html/2609.26972#bib.bib48),[29](https://arxiv.org/html/2609.26972#bib.bib49)\]\. These systems target classification and adapt final layers, and none trains a dynamics model whose gradient passes through an integrator\. Derivative\-free simultaneous perturbation stochastic approximation \(SPSA\)\[[33](https://arxiv.org/html/2609.26972#bib.bib12),[34](https://arxiv.org/html/2609.26972#bib.bib13)\]estimates a gradient from two randomly perturbed loss evaluations and therefore needs no backpropagation tape, but its gradient variance grows with the parameter count, limiting it to very small networks on microcontroller\-class hardware\. Jet matching is complementary and stronger where sensors stream densely: it retains exact analytic gradients of its loss with no zeroth\-order variance penalty, while remaining allocation\-free and solver\-free, and it scales to the hundreds of parameters where perturbation methods degrade\.

### 2\.3Derivative Estimation, System Identification, and Collocation

Fitting a vector field directly to estimated state derivatives,x˙≈Fθ​\(x\)\\dot\{x\}\\approx F\_\{\\theta\}\(x\), is the oldest idea in system identification\[[35](https://arxiv.org/html/2609.26972#bib.bib8),[36](https://arxiv.org/html/2609.26972#bib.bib56)\]\. Spline and generalized\-smoothing methods estimate the trajectory and its derivatives first and then regress the parameters\[[37](https://arxiv.org/html/2609.26972#bib.bib35),[38](https://arxiv.org/html/2609.26972#bib.bib36)\], symbolic regression distills closed\-form laws from derivative data\[[39](https://arxiv.org/html/2609.26972#bib.bib31)\], and sparse identification of nonlinear dynamics \(SINDy\) selects a few terms from a candidate library by sparse regression on numerically differentiated data\[[40](https://arxiv.org/html/2609.26972#bib.bib7),[41](https://arxiv.org/html/2609.26972#bib.bib37)\], with extensions to learned coordinates\[[42](https://arxiv.org/html/2609.26972#bib.bib29)\]and implicit dynamics\[[43](https://arxiv.org/html/2609.26972#bib.bib30)\]\. The known weakness of this entire family is what our baseline’sK=1K=1ablation exhibits: numerical differentiation amplifies measurement noise, and least squares against corrupted targets is biased\. Regularized differentiation methods\[[44](https://arxiv.org/html/2609.26972#bib.bib32),[45](https://arxiv.org/html/2609.26972#bib.bib33)\]and Savitzky–Golay filtering\[[7](https://arxiv.org/html/2609.26972#bib.bib6),[46](https://arxiv.org/html/2609.26972#bib.bib34)\]control the amplification through local polynomial fits, whose first\-derivative bias scales as𝒪⁡\(Δ​teff4\)\\mathcal\{O\}\(\\Delta t\_\{\\mathrm\{eff\}\}^\{4\}\)and noise as𝒪⁡\(σ/Δ​teff\)\\mathcal\{O\}\(\\sigma/\\Delta t\_\{\\mathrm\{eff\}\}\)at polynomial order four\.

Our work builds on this line in three ways:

1. 1\.We introduce a second\-order jet loss term with exact analytic gradients, using acceleration estimates to constrain the Jacobian of the vector field rather than only its value\.
2. 2\.We observe that the standard decimated SG pipeline discards\(m−1\)/m\(m\-1\)/mof the stream and correlates adjacent windows, and we recover the discarded independent samples with full\-rate phase sampling at negligible device cost\.
3. 3\.We introduce closed\-form on\-device rules that place the filter at its noise\-dependent optimum and reject windows where the polynomial model itself fails\. The quality gate is in the spirit of robust statistics\[[47](https://arxiv.org/html/2609.26972#bib.bib11)\], but acts on the targets’ diagnosable failure mode rather than on loss residuals, a distinction our double\-pendulum analysis shows to be decisive because a flexible network fits biased targets with small residuals\.

Unlike PINNs, which evaluate differential equations across continuous spatio\-temporal grids via automatic differentiation, jet matching processes discrete time\-series data using localized filtering and closed\-form gradient evaluations\. Table[1](https://arxiv.org/html/2609.26972#S2.T1)summarizes the key differences\.

Table 1:Taxonomy of differential equation learning frameworks\. Memory figures for the adjoint baseline reflect measured usage, while jet matching figures represent static memory allocations\.
### 2\.4Gradient Computation, Optimization, and Verification

The two gradient primitives of Section[3](https://arxiv.org/html/2609.26972#S3), the vector–Jacobian product and the Jacobian–vector product, are the elementary operations of reverse\- and forward\-mode automatic differentiation respectively\[[48](https://arxiv.org/html/2609.26972#bib.bib46),[49](https://arxiv.org/html/2609.26972#bib.bib43)\]\. Composing reverse mode over a forward\-mode pass is the classical technique for exact Hessian\-vector, and in this case Jacobian\-derivative products without materializing matrices\[[50](https://arxiv.org/html/2609.26972#bib.bib42)\]\. We verify both our gradients and the baseline’s adjoint against complex\-step differentiation, which evaluates a real function at a small imaginary perturbation and reads the derivative from the imaginary part, avoiding the subtractive cancellation of finite differences and attaining machine precision\[[51](https://arxiv.org/html/2609.26972#bib.bib44),[52](https://arxiv.org/html/2609.26972#bib.bib45)\]\. On the optimization side, our trainer combines stochastic approximation\[[53](https://arxiv.org/html/2609.26972#bib.bib40)\]in the form of Adam\[[9](https://arxiv.org/html/2609.26972#bib.bib9)\]with a cosine learning\-rate schedule\[[54](https://arxiv.org/html/2609.26972#bib.bib38)\]and Polyak–Ruppert iterate averaging\[[10](https://arxiv.org/html/2609.26972#bib.bib10)\], whose modern weight\-space form is known to find flatter minima and generalize better\[[55](https://arxiv.org/html/2609.26972#bib.bib39)\]\. The fixed\-memory sample buffer is maintained by reservoir sampling\[[56](https://arxiv.org/html/2609.26972#bib.bib41)\], which keeps a uniform random subset of an unbounded stream\. Each of these components is standard, and the contribution of this paper is showing that, arranged into a streaming, allocation\-free pipeline, they eliminate the accuracy penalty that has historically separated derivative\-matching from solver\-based training\.

## 3Methodology

### 3\.1Problem Formulation and Dynamical Systems

We consider a state observation stream given byx⁡\(tk\)\+ηk∈ℝnx\(t\_\{k\}\)\+\\eta\_\{k\}\\in\\mathbb\{R\}^\{n\}sampled at frequencyfs=1f\_\{s\}=1kHz, whereηk\\eta\_\{k\}represents independent and identically distributed \(i\.i\.d\.\) Gaussian sensor noise\. The true physical system followsx˙=f∗​\(x\)\\dot\{x\}=f^\{\\ast\}\(x\), which we decompose into:

f∗​\(x\)=fknown​\(x\)\+Δ⁡\(x\),f^\{\\ast\}\(x\)=f\_\{\\mathrm\{known\}\}\(x\)\+\\Delta\(x\),\(1\)wherefknownf\_\{\\mathrm\{known\}\}denotes the modeled physical dynamics andΔ\\Deltarepresents unknown dynamics\. The target learning model is structured as:

x˙=Fθ​\(x\)=fknown​\(x\)\+Nθ​\(x\),\\dot\{x\}=F\_\{\\theta\}\(x\)=f\_\{\\mathrm\{known\}\}\(x\)\+N\_\{\\theta\}\(x\),\(2\)whereNθN\_\{\\theta\}is a multilayer perceptron usingtanh\\tanhactivation functions\. The neural network architecture uses a2→16→16→22\{\\to\}16\{\\to\}16\{\\to\}2configuration \(354 parameters\) for the single pendulum and a4→16→16→44\{\\to\}16\{\\to\}16\{\\to\}4configuration \(420 parameters\) for the double pendulum\. Input and output layers are normalized using scale factors determined during an initial calibration phase\.

We evaluate performance on two physical systems:

1. 1\.Damped Single Pendulum:Defined by state vectorx=\[θ,ω\]⊤x=\[\\theta,\\omega\]^\{\\\!\\top\}with governing equations: θ˙=ω,\\dot\{\\theta\}=\\omega,\(3\)ω˙=−gL​sin⁡θ−c​ω,\\dot\{\\omega\}=\-\\frac\{g\}\{L\}\\sin\\theta\-c\\,\\omega,\(4\)whereg=9\.81​m/s2g=9\.81~\\mathrm\{m/s^\{2\}\},L=1​mL=1~\\mathrm\{m\}, and true damping coefficientc=0\.40c=0\.40\. The system is initialized from\(θ0,ω0\)=\(1\.2,0\)\(\\theta\_\{0\},\\omega\_\{0\}\)=\(1\.2,0\)with periodic random re\-excitation forces\.
2. 2\.Chaotic Damped Double Pendulum:Defined by state vectorx=\[θ1,θ2,ω1,ω2\]⊤x=\[\\theta\_\{1\},\\theta\_\{2\},\\omega\_\{1\},\\omega\_\{2\}\]^\{\\\!\\top\}, with link massesm1=m2=1m\_\{1\}=m\_\{2\}=1kg and lengthsl1=l2=1l\_\{1\}=l\_\{2\}=1m\. Defining the angular differenceδ=θ1−θ2\\delta=\\theta\_\{1\}\-\\theta\_\{2\}and common denominator termD=2​m1\+m2−m2​cos⁡\(2​δ\)\\mathrm\{D\}=2m\_\{1\}\+m\_\{2\}\-m\_\{2\}\\cos\(2\\delta\), the equations of motion are: θ˙1=ω1,\\dot\{\\theta\}\_\{1\}=\\omega\_\{1\},\(5\)θ˙2=ω2,\\dot\{\\theta\}\_\{2\}=\\omega\_\{2\},\(6\)ω˙1=−g\(2m1\+m2\)sinθ1−m2gsin\(θ1−2θ2\)−2sinδm2\(ω22l2\+ω12l1cosδ\)l1​D−c1​ω1,\\dot\{\\omega\}\_\{1\}=\\frac\{\-g\(2m\_\{1\}\{\+\}m\_\{2\}\)\\sin\\theta\_\{1\}\-m\_\{2\}g\\sin\(\\theta\_\{1\}\{\-\}2\\theta\_\{2\}\)\-2\\sin\\delta\\,m\_\{2\}\\big\(\\omega\_\{2\}^\{2\}l\_\{2\}\+\\omega\_\{1\}^\{2\}l\_\{1\}\\cos\\delta\\big\)\}\{l\_\{1\}\\,\\mathrm\{D\}\}\-c\_\{1\}\\,\\omega\_\{1\},\(7\)ω˙2=2sinδ\(ω12l1\(m1\+m2\)\+g\(m1\+m2\)cosθ1\+ω22l2m2cosδ\)l2​D−c2​ω2,\\dot\{\\omega\}\_\{2\}=\\frac\{2\\sin\\delta\\,\\big\(\\omega\_\{1\}^\{2\}l\_\{1\}\(m\_\{1\}\{\+\}m\_\{2\}\)\+g\(m\_\{1\}\{\+\}m\_\{2\}\)\\cos\\theta\_\{1\}\+\\omega\_\{2\}^\{2\}l\_\{2\}m\_\{2\}\\cos\\delta\\big\)\}\{l\_\{2\}\\,\\mathrm\{D\}\}\-c\_\{2\}\\,\\omega\_\{2\},\(8\)with damping parametersc1=c2=0\.150c\_\{1\}=c\_\{2\}=0\.150\. Initial conditions are set to\(2\.0,2\.4,0,0\)\(2\.0,2\.4,0,0\)within the chaotic regime\.

To evaluate model identification, we construct three UDE experimental configurations:

- •Case A:Damping terms are omitted fromfknownf\_\{\\mathrm\{known\}\}, requiringNθN\_\{\\theta\}to learn the full damping effect\.
- •Case B:Damping is included infknownf\_\{\\mathrm\{known\}\}, but gravity is mis\-specified asgmodel=0\.90​gg\_\{\\mathrm\{model\}\}=0\.90\\,g\. Here, the correct learned damping residual is zero\.
- •Case C:Both structural errors are present simultaneously \(mis\-specified gravity and omitted damping\)\.

Model accuracy is quantified using the relative vector field error over state domain𝒳\\mathcal\{X\}:

efield=∥Fθ−f∗∥L2​\(𝒳\)∥f∗∥L2​\(𝒳\),e\_\{\\mathrm\{field\}\}=\\frac\{\\lVert F\_\{\\theta\}\-f^\{\\ast\}\\rVert\_\{L\_\{2\}\(\\mathcal\{X\}\)\}\}\{\\lVert f^\{\\ast\}\\rVert\_\{L\_\{2\}\(\\mathcal\{X\}\)\}\},\(9\)evaluated over a41×4141\{\\times\}41grid for the single pendulum and over trajectories visited during evaluation for the double pendulum\. Recovered damping coefficientsc^\\hat\{c\}are computed by linear regression of the learned network outputs against state velocities\.

### 3\.2Lie–Taylor Jet Matching Loss and Analytic Gradients

For any smooth system following \([2](https://arxiv.org/html/2609.26972#S3.E2)\), the first and second time derivatives satisfyx˙=Fθ​\(x\)\\dot\{x\}=F\_\{\\theta\}\(x\)andx¨=JFθ​\(x\)​Fθ​\(x\)\\ddot\{x\}=J\_\{F\_\{\\theta\}\}\(x\)\\,F\_\{\\theta\}\(x\)via the chain rule\. These are the first terms of the Lie–Taylor expansion of the local flow, which is the Taylor series in time whose coefficients are successive Lie derivatives of the state along the vector field\[[8](https://arxiv.org/html/2609.26972#bib.bib55)\]\. Truncating it at second order gives the jet that the method matches\. Jet matching minimizes the discrepancy between these expressions and filtered derivative estimates\(x^,x˙^,x¨^\)\(\\hat\{x\},\\hat\{\\dot\{x\}\},\\hat\{\\ddot\{x\}\}\):

We define the first\-derivative residual vectorr1∈ℝnr\_\{1\}\\in\\mathbb\{R\}^\{n\}as:

r1=x˙^−Fθ​\(x^\)s1,r\_\{1\}=\\frac\{\\hat\{\\dot\{x\}\}\-F\_\{\\theta\}\(\\hat\{x\}\)\}\{s\_\{1\}\},\(10\)wheres1s\_\{1\}is the calibration Root Mean Square \(RMS\) scale for the first derivative\. Similarly, the second\-order model predictionGθ​\(x^\)∈ℝnG\_\{\\theta\}\(\\hat\{x\}\)\\in\\mathbb\{R\}^\{n\}is defined as:

Gθ​\(x^\)=JFθ​\(x^\)​Fθ​\(x^\),G\_\{\\theta\}\(\\hat\{x\}\)=J\_\{F\_\{\\theta\}\}\(\\hat\{x\}\)\\,F\_\{\\theta\}\(\\hat\{x\}\),\(11\)and the second\-derivative residual vectorr2∈ℝnr\_\{2\}\\in\\mathbb\{R\}^\{n\}is given by:

r2=x¨^−Gθ​\(x^\)s2,r\_\{2\}=\\frac\{\\hat\{\\ddot\{x\}\}\-G\_\{\\theta\}\(\\hat\{x\}\)\}\{s\_\{2\}\},\(12\)wheres2s\_\{2\}is the calibration RMS scale for the second derivative\. Combining these components yields the complete jet matching loss function:

ℒ⁡\(θ\)=λ1​∥r1∥22\+λ2​∥r2∥22,\\mathcal\{L\}\(\\theta\)=\\lambda\_\{1\}\\lVert r\_\{1\}\\rVert\_\{2\}^\{2\}\+\\lambda\_\{2\}\\lVert r\_\{2\}\\rVert\_\{2\}^\{2\},\(13\)whereλ1\\lambda\_\{1\}andλ2\\lambda\_\{2\}represent scalar loss weights\.

Because parametersθ\\thetaenter solely throughNθN\_\{\\theta\}, the exact per\-sample loss gradient can be computed analytically without numerical differential equation solvers\. We define the backpropagated cotangent vectoru¯∈ℝn\\bar\{u\}\\in\\mathbb\{R\}^\{n\}as:

u¯=−2​\(λ1​r1s1\+λ2​JFθ​\(x^\)⊤​r2s2\),\\bar\{u\}=\-2\\left\(\\lambda\_\{1\}\\,\\frac\{r\_\{1\}\}\{s\_\{1\}\}\+\\lambda\_\{2\}\\,J\_\{F\_\{\\theta\}\}\(\\hat\{x\}\)^\{\\\!\\top\}\\frac\{r\_\{2\}\}\{s\_\{2\}\}\\right\),\(14\)and the second\-order weight vectorw∈ℝnw\\in\\mathbb\{R\}^\{n\}as:

w=−2​r2s2\.w=\-2\\,\\frac\{r\_\{2\}\}\{s\_\{2\}\}\.\(15\)The full gradient of the loss function with respect toθ\\thetais then:

∇θℒ=\(∂Nθ​\(x^\)∂θ\)⊤​u¯\+λ2​∇θ⟨w,JNθ​\(x^\)​v⟩\|v=sg​\[Fθ​\(x^\)\],\\nabla\_\{\\theta\}\\mathcal\{L\}=\\left\(\\frac\{\\partial N\_\{\\theta\}\(\\hat\{x\}\)\}\{\\partial\\theta\}\\right\)^\{\\\!\\top\}\\bar\{u\}\+\\lambda\_\{2\}\\,\\nabla\_\{\\theta\}\\big\\langle w,\\;J\_\{N\_\{\\theta\}\}\(\\hat\{x\}\)\\,v\\big\\rangle\\Big\|\_\{\\,v=\\mathrm\{sg\}\[F\_\{\\theta\}\(\\hat\{x\}\)\]\},\(16\)wheresg⁡\[⋅\]\\mathrm\{sg\}\[\\cdot\]denotes a stop\-gradient operator holding vectorvvconstant during derivative evaluation\. This formulation needs only two primitive operations: a vector–Jacobian product \(VJP, a reverse\-mode pass that propagates an output cotangent back to the parameters\) through the network for the first term, and reverse\-mode differentiation over a Jacobian–vector product \(JVP, a forward\-mode pass that propagates a tangent directionvvthrough the network\) for the second term\[[49](https://arxiv.org/html/2609.26972#bib.bib43),[50](https://arxiv.org/html/2609.26972#bib.bib42)\]\. Both are evaluated within fixed layer buffers, so then×pn\\times pJacobian∂Nθ/∂θ\\partial N\_\{\\theta\}/\\partial\\thetais never formed and no dynamic memory is allocated during training\. Detailed mathematical derivations are provided in Appendix[A](https://arxiv.org/html/2609.26972#A1)\.

### 3\.3Streaming Derivative Estimation: Full\-Rate Phase Sampling

Derivative targets are derived from an1111\-tap \(M=5M=5\),4​th4\\mathrm\{th\}\-order polynomial \(P=4P=4\) Savitzky–Golay filter operating with tap stridemm\. The effective sampling step is given by:

Δ​teff=mfs\.\\Delta t\_\{\\mathrm\{eff\}\}=\\frac\{m\}\{f\_\{s\}\}\.\(17\)Using relative tap indicesτ∈\{−M,…,M\}\\tau\\in\\\{\-M,\\dots,M\\\}, we construct the Vandermonde matrixA∈ℝ\(2​M\+1\)×\(P\+1\)A\\in\\mathbb\{R\}^\{\(2M\+1\)\\times\(P\+1\)\}\(the design matrix of the local polynomial fit, whose columns are the monomialsτj\\tau^\{j\}evaluated at the tap positions\) with elements:

Ai​j=τij,i∈\{1,…,2​M\+1\},j∈\{0,…,P\}\.A\_\{ij\}=\\tau\_\{i\}^\{\\,j\},\\qquad i\\in\\\{1,\\dots,2M\+1\\\},\\quad j\\in\\\{0,\\dots,P\\\}\.\(18\)The row vectorsDdD\_\{d\}for derivative orderd∈\{0,1,2\}d\\in\\\{0,1,2\\\}are given by:

Dd=d\!\[\(A⊤A\)−1A⊤\]d:,D\_\{d\}=d\!\\,\\big\[\(A^\{\\\!\\top\}A\)^\{\-1\}A^\{\\\!\\top\}\\big\]\_\{d\\,:\},\(19\)where indexd:d\\,:selects the corresponding row of the pseudoinverse\. Given a sampled data vector𝐲∈ℝ2​M\+1\\mathbf\{y\}\\in\\mathbb\{R\}^\{2M\+1\}, state and derivative estimates are calculated as:

x^=D0⊤​𝐲,\\hat\{x\}=D\_\{0\}^\{\\\!\\top\}\\mathbf\{y\},\(20\)x˙^=D1⊤​𝐲Δ​teff,\\hat\{\\dot\{x\}\}=\\frac\{D\_\{1\}^\{\\\!\\top\}\\mathbf\{y\}\}\{\\Delta t\_\{\\mathrm\{eff\}\}\},\(21\)x¨^=D2⊤​𝐲Δ​teff2\.\\hat\{\\ddot\{x\}\}=\\frac\{D\_\{2\}^\{\\\!\\top\}\\mathbf\{y\}\}\{\\Delta t\_\{\\mathrm\{eff\}\}^\{2\}\}\.\(22\)
Standard decimation updates filtering windows everymmsteps, discarding intermediate samples and generating correlated outputs across adjacent windows\. To address this, we implement full\-rate phase sampling by sliding the filtering window by a small strides=5s=5while maintaining tap spacingmm\. Because interleaved window phases do not share raw measurement samples, their noise distributions are statistically independent\. This increases the effective volume of unique training targets by approximatelym/sm/swithout requiring additional sensors\. The required buffer footprint is limited to11​m11mfloating\-point values per channel \(4\.44\.4kB for single pendulum experiments atm=50m=50\)\. Figure[1](https://arxiv.org/html/2609.26972#S3.F1)illustrates this sampling approach\. Filtered outputs are collected into a reservoir buffer of sizeBB\(B=4096B=4096for single pendulum,B=2048B=2048for double pendulum\), from which minibatches of size6464are drawn for optimization\.

Figure 1:Overview of standard decimated sampling versus full\-rate phase sampling\. Maintaining tap spacingmmwhile sliding the filter window by stridessyieldsm/sm/sindependent derivative estimation streams\.
### 3\.4Noise\-Adaptive Filter Tuning and Quality Gating

We employ three adaptive mechanisms to automatically adjust filter settings based on operational noise conditions:

On\-Device Noise Estimation\.LetD0D\_\{0\}represent the smoothing row vector andμ\\mudenote the center tap index\. For a signal corrupted by zero\-mean Gaussian noise with varianceσ2\\sigma^\{2\}, the scalar residuale=yμ−D0⊤​𝐲e=y\_\{\\mu\}\-D\_\{0\}^\{\\\!\\top\}\\mathbf\{y\}has expected varianceσ2​\(1−2​\[D0\]μ\+∥D0∥22\)\\sigma^\{2\}\(1\-2\[D\_\{0\}\]\_\{\\mu\}\+\\lVert D\_\{0\}\\rVert\_\{2\}^\{2\}\)\. Averaging squared residuals overWWsliding windows yields an estimator for noise variance:

σ^2=1W​∑k=1Wek21−2​\[D0\]μ\+∥D0∥22,\\hat\{\\sigma\}^\{2\}=\\frac\{\\frac\{1\}\{W\}\\sum\_\{k=1\}^\{W\}e\_\{k\}^\{\\,2\}\}\{1\-2\\,\[D\_\{0\}\]\_\{\\mu\}\+\\lVert D\_\{0\}\\rVert\_\{2\}^\{2\}\},\(23\)which incurs minimal overhead \(one multiply\-accumulate operation per tap\)\.

Adaptive Time\-Scale Selection\.Balancing polynomial approximation bias𝒪⁡\(Δ​teff4\)\\mathcal\{O\}\(\\Delta t\_\{\\mathrm\{eff\}\}^\{4\}\)against noise variance𝒪⁡\(σ/Δ​teff\)\\mathcal\{O\}\(\\sigma/\\Delta t\_\{\\mathrm\{eff\}\}\)yields an optimal time stepΔ​teff∗∝σ1/5\\Delta t\_\{\\mathrm\{eff\}\}^\{\\ast\}\\propto\\sigma^\{1/5\}\. Scaling this relationship inversely with the characteristic signal frequencyωc=s2/s1\\omega\_\{c\}=s\_\{2\}/s\_\{1\}gives:

Δ​teff∗=0\.09ωc​\(σ^rel0\.01\)1/5,\\Delta t\_\{\\mathrm\{eff\}\}^\{\\ast\}=\\frac\{0\.09\}\{\\omega\_\{c\}\}\\left\(\\frac\{\\hat\{\\sigma\}\_\{\\mathrm\{rel\}\}\}\{0\.01\}\\right\)^\{1/5\},\(24\)whereσ^rel\\hat\{\\sigma\}\_\{\\mathrm\{rel\}\}is the estimated relative noise standard deviation\. The constant0\.090\.09serves as a baseline scaling factor determined from single pendulum calibration experiments\.

Second\-Order Derivative Gating\.Second derivative estimatesx¨^\\hat\{\\ddot\{x\}\}carry higher noise standard deviations \(∥D2∥2​σ^/Δ​teff2\\lVert D\_\{2\}\\rVert\_\{2\}\\hat\{\\sigma\}/\\Delta t\_\{\\mathrm\{eff\}\}^\{2\}\)\. To prevent noise amplification from degrading optimization, we compute the channel Signal\-to\-Noise Ratio \(SNR\):

SNR2=minchannels⁡s2∥D2∥2​σ^/Δ​teff2\.\\mathrm\{SNR\}\_\{2\}=\\min\_\{\\text\{channels\}\}\\frac\{s\_\{2\}\}\{\\lVert D\_\{2\}\\rVert\_\{2\}\\,\\hat\{\\sigma\}/\\Delta t\_\{\\mathrm\{eff\}\}^\{2\}\}\.\(25\)The second\-order loss weightλ2\\lambda\_\{2\}is adaptively assigned based on thresholdτSNR=12\\tau\_\{\\mathrm\{SNR\}\}=12:

λ2=\{λ20if​SNR2≥τSNR,0if​SNR2<τSNR,\\lambda\_\{2\}=\\begin\{cases\}\\lambda\_\{2\}^\{\\,0\}&\\text\{if \}\\mathrm\{SNR\}\_\{2\}\\geq\\tau\_\{\\mathrm\{SNR\}\},\\\\\[2\.0pt\] 0&\\text\{if \}\\mathrm\{SNR\}\_\{2\}<\\tau\_\{\\mathrm\{SNR\}\},\\end\{cases\}\(26\)whereλ20=0\.1\\lambda\_\{2\}^\{\\,0\}=0\.1\. IfSNR2\\mathrm\{SNR\}\_\{2\}drops below the threshold, second\-order derivative matching is disabled to protect gradient updates from high variance\.

Polynomial Misfit Quality Gating\.High\-frequency non\-linear dynamics can violate local polynomial assumptions during rapid state transitions\. To detect these conditions, we form the residual projection matrixR=I−A​A\+R=I\-AA^\{\+\}, whereA\+A^\{\+\}is the pseudoinverse ofAA\. The fit residual norm∥R​𝐲∥2\\lVert R\\,\\mathbf\{y\}\\rVert\_\{2\}quantifies polynomial deviation\. Samples are rejected if they exceed a dynamic threshold:

∥R​𝐲∥2\>6⋅median⁡\(∥R​𝐲∥2\)\+0\.02⋅RMS⁡\(𝐲\),\\lVert R\\,\\mathbf\{y\}\\rVert\_\{2\}\>6\\cdot\\mathrm\{median\}\\big\(\\lVert R\\,\\mathbf\{y\}\\rVert\_\{2\}\\big\)\+0\.02\\cdot\\mathrm\{RMS\}\(\\mathbf\{y\}\),\(27\)where the median is computed over recent windows and the RMS term establishes a baseline floor\. Rejected samples are excluded from the training buffer\. Algorithm[1](https://arxiv.org/html/2609.26972#alg1)outlines the complete training procedure\.

Algorithm 1Noise\-Adaptive Lie–Taylor Jet Matching Routine1:Sensor stream

yky\_\{k\}at

fs=1f\_\{s\}=1kHz; baseline model

fknownf\_\{\\mathrm\{known\}\}; parameters from Table[5](https://arxiv.org/html/2609.26972#A2.T5)

2:Optimized weights

θEMA\\theta\_\{\\mathrm\{EMA\}\}
3:Calibration:Measure scales

xrms,s1,s2x\_\{\\mathrm\{rms\}\},s\_\{1\},s\_\{2\}using SG estimates \([20](https://arxiv.org/html/2609.26972#S3.E20)\)–\([22](https://arxiv.org/html/2609.26972#S3.E22)\)

4:Adaptation:Estimate noise

σ^\\hat\{\\sigma\}using \([23](https://arxiv.org/html/2609.26972#S3.E23)\); set

ωc←s2/s1\\omega\_\{c\}\\leftarrow s\_\{2\}/s\_\{1\},

m←fs​Δ​teff∗m\\leftarrow f\_\{s\}\\Delta t\_\{\\mathrm\{eff\}\}^\{\\ast\}via \([24](https://arxiv.org/html/2609.26972#S3.E24)\), and weight

λ2\\lambda\_\{2\}via \([26](https://arxiv.org/html/2609.26972#S3.E26)\)

5:foreach incoming sample tick

kkdo

6:Store

yky\_\{k\}in circular buffer of length

\(2​M\+1\)​m\(2M\{\+\}1\)mper channel

7:if

kmods==0k\\bmod s==0then

8:Construct window vector

𝐲\\mathbf\{y\}; evaluate estimates

\(x^,x˙^,x¨^\)\(\\hat\{x\},\\hat\{\\dot\{x\}\},\\hat\{\\ddot\{x\}\}\)via \([20](https://arxiv.org/html/2609.26972#S3.E20)\)–\([22](https://arxiv.org/html/2609.26972#S3.E22)\)

9:if

∥R​𝐲∥2\\lVert R\\,\\mathbf\{y\}\\rVert\_\{2\}satisfies condition \([27](https://arxiv.org/html/2609.26972#S3.E27)\) across all channelsthen

10:Add sample

\(x^,x˙^,x¨^\)\(\\hat\{x\},\\hat\{\\dot\{x\}\},\\hat\{\\ddot\{x\}\}\)to reservoir buffer of capacity

BB
11:endif

12:endif

13:if

kmod\(fs/ftrain\)==0k\\bmod\(f\_\{s\}/f\_\{\\mathrm\{train\}\}\)==0andbuffer is filledthen

14:Draw minibatch of 64 samples; evaluate gradient

∇θℒ\\nabla\_\{\\theta\}\\mathcal\{L\}using \([14](https://arxiv.org/html/2609.26972#S3.E14)\)–\([16](https://arxiv.org/html/2609.26972#S3.E16)\)

15:Apply Adam update with cosine learning rate schedule; clip gradient norm at

5\.05\.0
16:Update exponential moving average:

θEMA←β​θEMA\+\(1−β\)​θ\\theta\_\{\\mathrm\{EMA\}\}\\leftarrow\\beta\\theta\_\{\\mathrm\{EMA\}\}\+\(1\{\-\}\\beta\)\\theta
17:endif

18:endfor

19:Continue training updates on stored reservoir buffer until reaching total iteration limit

20:return

θEMA\\theta\_\{\\mathrm\{EMA\}\}

### 3\.5Static Memory Allocation

The training algorithm relies entirely on static memory allocations: network parameters, gradient buffers, Adam optimizer states, and EMA copies require5×3545\\times 354floating\-point values for the single pendulum\. Including the sample reservoir \(B×3​nB\\times 3n\), circular input buffers \(11​m11mvalues/channel\), filter coefficients, scale vectors, and projection matrixRR, total static memory usage is108\.1108\.1kB\. This memory footprint remains constant throughout execution, avoiding dynamic allocations during online training\. Section[5\.6](https://arxiv.org/html/2609.26972#S5.SS6)reports an on\-device measurement of this footprint at a reduced reservoir capacity \(B=2048B=2048, the same capacity used for the double pendulum in Table[5](https://arxiv.org/html/2609.26972#A2.T5)\), adopted to fit the target board’s available DRAM alongside the Arduino/ESP\-IDF runtime\.

## 4Experimental Setup

### 4\.1Adjoint Baseline Configuration

We compare the proposed method against a standard solver\-based UDE training pipeline: fixed\-step 4th\-order Runge–Kutta \(RK4\) integration inside the loss function, multiple shooting overKK\-step windows initialized from trajectory data\[[3](https://arxiv.org/html/2609.26972#bib.bib5)\], exact hand\-derived discrete adjoint backpropagation across integration steps \(reverse\-mode differentiation applied to the discretized solver steps, as opposed to the continuous adjoint\), and full\-batch Adam optimization\[[9](https://arxiv.org/html/2609.26972#bib.bib9)\]\. The discrete adjoint implementation was validated against complex\-step differentiation, matching numerical gradients to within1\.8×10−151\.8\{\\times\}10^\{\-15\}to5\.2×10−135\.2\{\\times\}10^\{\-13\}across test cases and shooting window lengthsK∈\{1,3,5,10\}K\\in\\\{1,3,5,10\\\}\. The baseline reference configuration \(K=40K=40, 128 windows, 1,200 iterations, 50 Hz observation rate\) requires a5\.95\.9MB intermediate tape,9\.19\.1MB peak heap usage, and26\.226\.2million multiply\-accumulate \(MAC\) operations per iteration\.

### 4\.2Experimental Control and Data Matching

Both jet matching and adjoint baseline implementations share common evaluation components: physical dynamics routines, synthetic noise generation, initial conditions, network architectures, weight initializations, and field error calculation code\. Experimental conditions are matched across tests using identical dataset durations \(6060s for single pendulum,2020s for double pendulum\)\. Sensitivity checks were performed on baseline parameters: increasing the baseline observation sampling rate to 100 Hz resulted in higher error \(0\.00330\.0033vs\.0\.00280\.0028at the reference5050Hz rate, 1% noise\), while increasing the shooting window toK=80K=80doubled tape size to11\.811\.8MB\. Single pendulum experiments were evaluated across 3 random seeds, and double pendulum experiments across 4 seeds\.

### 4\.3Execution Environment and Resource Projections

Accuracy evaluations were conducted using 64\-bit floating\-point simulations in Python/NumPy to isolate algorithmic behavior under controlled numerical precision\. Microcontroller hardware projections target an ESP32 microcontroller \(Xtensa LX6 dual\-core processor running at 240 MHz, 320 kB internal SRAM, single\-precision floating\-point unit \(FPU\)\)\. Memory figures reflect static compiler allocations\. Computation times are projected based on counted operations \(184,960184\{,\}960MACs per update cycle\) assuming an execution throughput of 3 clock cycles per MAC and 60 cycles per transcendental function call \(tanhf\), compared against desktop baseline measurements\. Section[5\.6](https://arxiv.org/html/2609.26972#S5.SS6)replaces this projection with a direct on\-device measurement from the same training algorithm running on physical ESP32 hardware\.

## 5Experiments and Results

### 5\.1Accuracy Across Sensor Noise Levels

Table[2](https://arxiv.org/html/2609.26972#S5.T2)presents relative field error results across varying sensor noise levels\. An unadapted jet matching baseline \(v1: decimated sampling, minibatch size 32, fixed learning rates, constantλ2\\lambda\_\{2\}\) produces a relative field error of0\.00900\.0090at 1% noise, a3\.2×3\.2\\timeshigher error than the adjoint baseline \(0\.00280\.0028\)\. Incorporating the adaptive mechanisms described in Section[3](https://arxiv.org/html/2609.26972#S3)\(v2\) reduces the relative field error to0\.00220\.0022under identical 60 s data windows\. Across noise levels from 0% to 5%, the proposed method achieves lower mean field error than the reference baseline in five out of six evaluated conditions, yielding a geometric\-mean error ratio of0\.65×0\.65\\times\(0\.82×0\.82\\timescompared to the baseline usingK=80K=80\)\. Parameter estimation accuracy remains consistent across noise levels, recovering single pendulum damping parameters withinc^∈\[0\.388,0\.399\]\\hat\{c\}\\in\[0\.388,0\.399\]\(true value0\.400\.40\), while baseline estimates drift to0\.4340\.434under 5% noise\.

Table 2:Relative field error \(μ±σ\\mu\\pm\\sigmaover 3 seeds\), single pendulum Case A, matched 60 s data window\. Bold text indicates lowest error per noise level\.†\\dagger: 1–3 seeds\. Error ratios<1\.0<1\.0indicate lower relative error for jet matching\.Figure[2](https://arxiv.org/html/2609.26972#S5.F2)illustrates the solver\-based training loop whose solver and tape jet matching removes\. Table[3](https://arxiv.org/html/2609.26972#S5.T3)summarizes the operational head\-to\-head comparison, while Figure[3](https://arxiv.org/html/2609.26972#S5.F3)illustrates the relationship between shooting window lengthKK, memory usage, and noise performance in the baseline model\.

Figure 2:The solver\-based training loop that jet matching removes\. Each gradient requires4​K4Kvector\-field evaluations forward and again in reverse, and the checkpoint tape retains every solver stage and activation for allKKsteps\. Jet matching replaces the entire dashed optimization loop and the tape with the streaming pipeline of Algorithm[1](https://arxiv.org/html/2609.26972#alg1)\.Table 3:Detailed comparison under matched physical plants, neural network structures, noise models, and evaluation routines\.\*Projected for theB=4096B=4096reference configuration \(Table[5](https://arxiv.org/html/2609.26972#A2.T5)\)\. Section[5\.6](https://arxiv.org/html/2609.26972#S5.SS6)measures the same algorithm on physical ESP32 hardware at a reduced reservoir \(B=2048B=2048, fit to the target board’s DRAM\):61\.361\.3kB static memory and7\.247\.24ms per update\.

![Refer to caption](https://arxiv.org/html/2609.26972v1/fig3_noise_shooting.png)Figure 3:Evaluation of the adjoint baseline: field error versus noise level atK=40K=40\(A\), accuracy across various shooting window lengthsKK\(B\), and trade\-offs between computational cost and field accuracy \(C\)\.
### 5\.2Evaluation on Chaotic Double Pendulum Dynamics

Under 1% sensor noise over a 20 s dataset, jet matching achieves a mean relative field error of0\.01150\.0115\(median0\.00830\.0083across 4 seeds\) on the chaotic double pendulum, compared to0\.01370\.0137for the adjoint baseline \(spread0\.00880\.0088to0\.01980\.0198across 3 seeds\)\. Evaluating jet matching models on trajectory states generated by the baseline yields a consistent field error of0\.01160\.0116\(median0\.00780\.0078\)\. The model recovers unmodeled damping coefficients atc^1=0\.141\\hat\{c\}\_\{1\}=0\.141\(true value0\.1500\.150\), compared to the baseline’s estimate of0\.1600\.160\.

Occasional error outliers were observed on this system, where 1 out of 4 seeds produced an error of0\.0250\.025\. Analysis indicates that during rapid chaotic flips, higher\-order Taylor series terms exceed the approximation capacity of 4th\-order polynomial filters over fixed window lengths\. This causes local derivative estimation errors that pass through the quality gate, temporarily distorting the learned field in high\-velocity states \(\|ω\|≫0\\lvert\\omega\\rvert\\gg 0\)\.

### 5\.3System Identification and Long\-Term Rollout Dynamics

Figure[4](https://arxiv.org/html/2609.26972#S5.F4)illustrates identification performance across different UDE structural formulations under 1% noise\. Case A recovers a damping coefficient ofc^=0\.397\\hat\{c\}=0\.397, Case C yieldsc^=0\.387\\hat\{c\}=0\.387, and Case B, where true damping is fully included infknownf\_\{\\mathrm\{known\}\}, returns an estimated residual coefficient of−0\.001\-0\.001, confirming that the network does not introduce unnecessary corrections when physical models are accurate\. Figure[5](https://arxiv.org/html/2609.26972#S5.F5)evaluates long\-term predictive stability\. When integrated over an88s trajectory window, the trained hybrid field maintains a trajectory error of0\.00680\.0068, representing a393×393\\timesreduction in error compared to unaugmented baseline models \(fknownf\_\{\\mathrm\{known\}\}alone\) and a10×10\\timesimprovement over unadapted jet matching implementations \(0\.0680\.068\)\.

Figure 4:System identification across UDE configurations under 1% noise: reconstructed residual functions versus true dynamics \(A, B\) and estimated damping parameters across random seeds \(C\)\.Figure 5:Autonomous trajectory rollouts over an 8 s horizon \(A\), phase space trajectories across multiple initial conditions \(B\), and relative state error progression over time \(C\)\.
### 5\.4Ablation Study and Performance Drivers

Table[4](https://arxiv.org/html/2609.26972#S5.T4)breaks down performance improvements by adding individual adaptive mechanisms under 1% noise\. The largest accuracy gains stem from full\-rate phase sampling and larger sample buffer sizes, which increase sample independence\. Dynamic time\-scale adaptation and quality gating provide additional refinements, particularly in low\-noise settings where fixed filter parameters introduce truncation errors\.

Table 4:Ablation study showing cumulative relative field error improvements on the single pendulum under 1% noise across 3 random seeds\.Figure[6](https://arxiv.org/html/2609.26972#S5.F6)illustrates the impact of second\-order derivative gating \(λ2\\lambda\_\{2\}\)\. Under 1% noise, the second\-order loss component remains pinned near its noise floor \(∼2×10−2\\sim 2\{\\times\}10^\{\-2\}\) for the entire run, while the first\-order loss component drops to1\.8×10−31\.8\{\\times\}10^\{\-3\}\. Disabling second\-order matching when SNR is low prevents high\-variance derivative estimates from impacting parameter updates\. Figure[7](https://arxiv.org/html/2609.26972#S5.F7)shows performance ratios across different noise levels and filter settings, confirming that second\-order matching provides minimal benefit when first\-order sample density is sufficient\.

Figure 6:Evaluation of derivative loss components during training at 1% noise\. The second\-order loss term remains near a noise floor \(A\), and disabling it yields equivalent performance \(B\), supporting the SNR gating strategy\.![Refer to caption](https://arxiv.org/html/2609.26972v1/fig7_noise_lam2.png)Figure 7:Performance comparison of second\-order loss contributions across noise levels and filter decimation steps: field error ratios \(λ2=0\.1\\lambda\_\{2\}=0\.1vs\.λ2=0\\lambda\_\{2\}=0\) relative to calculated SNR thresholds\.
### 5\.5Resource Usage and Memory Requirements

The baseline model requires5\.95\.9MB of intermediate tape storage \(11\.811\.8MB atK=80K=80\), peaks at9\.19\.1MB of heap allocation, and processes26\.226\.2M MACs per iteration\. In contrast, jet matching operates within a static memory allocation of108\.1108\.1kB, representing a5555to109×109\\timesreduction in memory usage\. Per\-update computational requirements are184,960184\{,\}960MACs \(batch size 64\), corresponding to a cycle\-counted estimate of2\.852\.85ms per step on a 240 MHz ESP32 microcontroller \(7\.1%7\.1\\%duty cycle at 25 Hz\)\. Figure[8](https://arxiv.org/html/2609.26972#S5.F8)C summarizes this projection alongside the corresponding measurement from physical ESP32 hardware, detailed in Section[5\.6](https://arxiv.org/html/2609.26972#S5.SS6)\.

Figure 8:Ablation analysis under 1% noise \(A\), noise\-free parameter recovery performance \(B\), and the ESP32 resource account, projected alongside the physical measurement of Section[5\.6](https://arxiv.org/html/2609.26972#S5.SS6)\(C\)\.
### 5\.6On\-Device Validation on ESP32 Hardware

We implemented the complete training algorithm, including all mechanisms of Section[3](https://arxiv.org/html/2609.26972#S3), in single\-precision C on an ESP32 \(dual\-core Xtensa LX6, 240 MHz, no floating\-point coprocessor for transcendental functions\) and ran it end\-to\-end on physical hardware at1%1\\%sensor noise on the single pendulum\. The target board’s available DRAM could not accommodate the reference reservoir capacity \(B=4096B=4096, Table[5](https://arxiv.org/html/2609.26972#A2.T5)\) alongside the Arduino/ESP\-IDF runtime’s own static allocations, so the reservoir was set toB=2048B=2048, the identical capacity already used for the double pendulum in simulation\. The on\-device noise estimator measuredσ^rel=0\.0086\\hat\{\\sigma\}\_\{\\mathrm\{rel\}\}=0\.0086and selected a decimation ofm=26m=26\(versusm≈39m\\approx 39in the float64 simulation at the same noise level\)\. Equation \([24](https://arxiv.org/html/2609.26972#S3.E24)\) depends only onωc\\omega\_\{c\}andσ^rel\\hat\{\\sigma\}\_\{\\mathrm\{rel\}\}, not on reservoir capacity, so this shift reflects estimation variance inωc\\omega\_\{c\}andσ^rel\\hat\{\\sigma\}\_\{\\mathrm\{rel\}\}over the shorter on\-device calibration window together with single\-precision arithmetic, rather than the change inBBitself\. The second\-order term remained gated off \(SNR2=2\.5<τSNR\\mathrm\{SNR\}\_\{2\}=2\.5<\\tau\_\{\\mathrm\{SNR\}\}\), consistent with Section[4](https://arxiv.org/html/2609.26972#S4)\.

Figure[9](https://arxiv.org/html/2609.26972#S5.F9)reports the result\. The on\-device run reaches a relative field error of0\.00200\.0020and recovers the damping coefficient toc^=0\.400\\hat\{c\}=0\.400\(true0\.4000\.400\), matching the accuracy of the float64 simulation at the corresponding operating point despite the smaller reservoir and reduced numerical precision; the quality gate rejected6\.5%6\.5\\%of candidate windows \(780780of12,00012\{,\}000\), consistent with the rate reported in Appendix[C](https://arxiv.org/html/2609.26972#A3)\. Measured static memory usage is61\.361\.3kB, exactly matching the value computed fromsizeofon the deployed arrays atB=2048B=2048, and a96×96\\timesreduction relative to the adjoint baseline’s5\.95\.9MB tape \(consistent with the5555–109×109\\timesrange of Section[5](https://arxiv.org/html/2609.26972#S5)\)\.

The measured update time is7\.247\.24ms, against the2\.852\.85ms cycle\-counted projection of Section[4](https://arxiv.org/html/2609.26972#S4), a gap we attribute to the projection’s idealized cost model: it counts multiply\-accumulates andtanhfcalls but not the floating\-point divisions the algorithm performs per sample \(normalizing bysxs\_\{x\},s1s\_\{1\},s2s\_\{2\}\) or ordinary loop and function\-call overhead the compiler does not fully eliminate\. The measured time is highly consistent across all8,0008\{,\}000updates of the run \(interquartile range under0\.20\.2ms; inset of Figure[9](https://arxiv.org/html/2609.26972#S5.F9)B\), indicating that the gap is a fixed constant\-factor cost rather than variable jitter\. At7\.247\.24ms per update the training loop occupies an18\.1%18\.1\\%duty cycle at the paper’s2525Hz training rate, still comfortably within budget\. We report the measured value as the operative timing figure and retain the projection only for its decomposition into MAC and transcendental\-call counts\.

Figure 9:On\-device validation on ESP32 hardware at1%1\\%noise,B=2048B=2048: relative field error and recoveredc^\\hat\{c\}over training \(A\), measured per\-update latency against the projected value and the2525Hz timing budget, with a zoomed inset on the measured spread \(B\), and the measured resource account against the simulated projection \(C\)\.

## 6Conclusion

### 6\.1Summary

This work evaluates whether numerical differential equation solvers can be removed from UDE training pipelines without sacrificing model accuracy\. By replacing internal solver integrations with noise\-adaptive Lie–Taylor jet matching, we achieve comparable or superior vector field accuracy across multiple test systems while reducing training memory requirements by roughly two orders of magnitude \(5555to109×109\\times\)\. Full\-rate phase sampling, reservoir buffering, adaptive filter tuning, and quality gating allow direct derivative matching to maintain robust performance under measurement noise\. The resulting approach operates within108108kB of statically allocated memory, enabling online physics\-informed model learning on edge microcontrollers\.

### 6\.2Limitations and Future Work

We highlight several operational constraints of the current approach:

1. 1\.We validated the complete algorithm on physical ESP32 hardware \(Section[5\.6](https://arxiv.org/html/2609.26972#S5.SS6)\) at a single operating point \(1%1\\%noise, single pendulum,B=2048B=2048\)\. Extending on\-device measurement across the full noise sweep, to the double pendulum, and to power consumption remains future work\.
2. 2\.In low\-noise scenarios \(0\.1%0\.1\\%noise\), traditional adjoint methods maintain an accuracy advantage \(1\.8×1\.8\\timeslower absolute error\)\.
3. 3\.On complex chaotic systems, rapid state transitions can occasionally exceed polynomial filter assumptions, requiring further refinements such as variable window lengths or dynamic filter ordering\.
4. 4\.Testing has focused on simulated physical systems under Gaussian noise assumptions, so evaluating performance on hardware platforms with non\-Gaussian noise, sensor latency, and quantization effects remains future work\.

## 7Acknowledgments

The authors acknowledge the use of AI tools during the preparation of this manuscript for assistance with rephrasing and paraphrasing text\. All AI\-generated content was thoroughly reviewed, verified, and edited by the authors, who take full responsibility for the accuracy and integrity of the work presented\.

## References

- \[1\]C\. Rackauckas, Y\. Ma, J\. Martensen, C\. Warner, K\. Zubov, R\. Supekar, D\. Skinner, A\. Ramadhan, and A\. Edelman\(2021\)Universal differential equations for scientific machine learning\.Note:arXivExternal Links:2001\.04385v4,[Link](https://arxiv.org/abs/2001.04385v4)Cited by:[§1\.1](https://arxiv.org/html/2609.26972#S1.SS1.p1.1),[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p1.1),[Table 1](https://arxiv.org/html/2609.26972#S2.T1.2.2.1)\.
- \[2\]\(2018\)Neural ordinary differential equations\.InAdvances in Neural Information Processing Systems \(NeurIPS\),Vol\.31,pp\. 6571–6583\.External Links:[Link](https://proceedings.neurips.cc/paper/2018/hash/69386f6bb1dfed68692a24c8686939b9-Abstract.html)Cited by:[§1\.1](https://arxiv.org/html/2609.26972#S1.SS1.p2.1),[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p1.1)\.
- \[3\]H\. G\. Bock and K\. J\. Plitt\(1984\)A multiple shooting algorithm for direct solution of optimal control problems\.IFAC Proceedings Volumes17\(2\),pp\. 1603–1608\.Note:9th IFAC World Congress: A Bridge Between Control Science and Technology, Budapest, Hungary, 2–6 July 1984External Links:ISSN 1474\-6670,[Document](https://dx.doi.org/10.1016/S1474-6670%2817%2961205-9),[Link](https://www.sciencedirect.com/science/article/pii/S1474667017612059)Cited by:[§1\.1](https://arxiv.org/html/2609.26972#S1.SS1.p2.1),[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p2.1),[§4\.1](https://arxiv.org/html/2609.26972#S4.SS1.p1.1)\.
- \[4\]R\. David, J\. Duke, A\. Jain, V\. J\. Reddi, N\. Jeffries, J\. Li, N\. Kreeger, I\. Nappier, M\. Natraj, S\. Regev, R\. Rhodes, T\. Wang, and P\. Warden\(2021\)TensorFlow Lite Micro: embedded machine learning on TinyML systems\.Note:arXivExternal Links:2010\.08678v3,[Link](https://arxiv.org/abs/2010.08678v3)Cited by:[item 1](https://arxiv.org/html/2609.26972#S1.I1.i1.p1.1),[§2\.2](https://arxiv.org/html/2609.26972#S2.SS2.p1.1)\.
- \[5\]J\. Lin, W\. Chen, Y\. Lin, J\. Cohn, C\. Gan, and S\. Han\(2020\)MCUNet: tiny deep learning on IoT devices\.Note:arXivExternal Links:2007\.10319,[Link](https://arxiv.org/abs/2007.10319)Cited by:[item 1](https://arxiv.org/html/2609.26972#S1.I1.i1.p1.1),[§2\.2](https://arxiv.org/html/2609.26972#S2.SS2.p1.1)\.
- \[6\]J\. Lin, L\. Zhu, W\. Chen, W\. Wang, C\. Gan, and S\. Han\(2024\)On\-device training under 256KB memory\.Note:arXivExternal Links:2206\.15472v4,[Link](https://arxiv.org/abs/2206.15472v4)Cited by:[item 1](https://arxiv.org/html/2609.26972#S1.I1.i1.p1.1),[§2\.2](https://arxiv.org/html/2609.26972#S2.SS2.p1.1)\.
- \[7\]J\. Steinier, Y\. Termonia, and J\. Deltour\(1972\)Smoothing and differentiation of data by simplified least square procedure\.Analytical Chemistry44\(11\),pp\. 1906–1909\.External Links:ISSN 0003\-2700,[Document](https://dx.doi.org/10.1021/ac60319a045),[Link](https://doi.org/10.1021/ac60319a045)Cited by:[1st item](https://arxiv.org/html/2609.26972#S1.I2.i1.p1.1),[§2\.3](https://arxiv.org/html/2609.26972#S2.SS3.p1.1)\.
- \[8\]A\. Isidori\(1995\)Nonlinear control systems\.3 edition,Communications and Control Engineering,Springer London\.External Links:ISBN 978\-3\-540\-19916\-8,[Document](https://dx.doi.org/10.1007/978-1-84628-615-5)Cited by:[1st item](https://arxiv.org/html/2609.26972#S1.I2.i1.p1.1),[§3\.2](https://arxiv.org/html/2609.26972#S3.SS2.p1.1)\.
- \[9\]D\. P\. Kingma and J\. Ba\(2015\)Adam: a method for stochastic optimization\.In3rd International Conference on Learning Representations \(ICLR\),External Links:[Document](https://dx.doi.org/10.48550/arXiv.1412.6980)Cited by:[item 3](https://arxiv.org/html/2609.26972#S1.I2.i2.I1.i3.p1.1),[§2\.4](https://arxiv.org/html/2609.26972#S2.SS4.p1.1),[§4\.1](https://arxiv.org/html/2609.26972#S4.SS1.p1.1)\.
- \[10\]B\. T\. Polyak and A\. B\. Juditsky\(1992\)Acceleration of stochastic approximation by averaging\.SIAM Journal on Control and Optimization30\(4\),pp\. 838–855\.External Links:[Document](https://dx.doi.org/10.1137/0330046)Cited by:[item 3](https://arxiv.org/html/2609.26972#S1.I2.i2.I1.i3.p1.1),[§2\.4](https://arxiv.org/html/2609.26972#S2.SS4.p1.1)\.
- \[11\]Y\. Rubanova, R\. T\. Q\. Chen, and D\. Duvenaud\(2019\)Latent ODEs for irregularly\-sampled time series\.Note:arXivExternal Links:1907\.03907,[Link](https://arxiv.org/abs/1907.03907)Cited by:[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p1.1)\.
- \[12\]E\. Dupont, A\. Doucet, and Y\. W\. Teh\(2019\)Augmented neural ODEs\.Note:arXivExternal Links:1904\.01681,[Link](https://arxiv.org/abs/1904.01681)Cited by:[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p1.1)\.
- \[13\]P\. Kidger\(2022\)On neural differential equations\.Note:arXivExternal Links:2202\.02435,[Link](https://arxiv.org/abs/2202.02435)Cited by:[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p1.1)\.
- \[14\]S\. Greydanus, M\. Dzamba, and J\. Yosinski\(2019\)Hamiltonian neural networks\.Note:arXivExternal Links:1906\.01563,[Link](https://arxiv.org/abs/1906.01563)Cited by:[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p1.1)\.
- \[15\]M\. Cranmer, S\. Greydanus, S\. Hoyer, P\. Battaglia, D\. Spergel, and S\. Ho\(2020\)Lagrangian neural networks\.Note:arXivExternal Links:2003\.04630,[Link](https://arxiv.org/abs/2003.04630)Cited by:[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p1.1)\.
- \[16\]M\. Lutter, C\. Ritter, and J\. Peters\(2019\)Deep Lagrangian networks: using physics as model prior for deep learning\.Note:arXivExternal Links:1907\.04490,[Link](https://arxiv.org/abs/1907.04490)Cited by:[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p1.1)\.
- \[17\]M\. Raissi, P\. Perdikaris, and G\. E\. Karniadakis\(2019\)Physics\-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations\.Journal of Computational Physics378,pp\. 686–707\.External Links:[Document](https://dx.doi.org/10.1016/j.jcp.2018.10.045)Cited by:[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p1.1),[Table 1](https://arxiv.org/html/2609.26972#S2.T1.2.3.1)\.
- \[18\]G\. E\. Karniadakis, I\. G\. Kevrekidis, L\. Lu, P\. Perdikaris, S\. Wang, and L\. Yang\(2021\)Physics\-informed machine learning\.Nature Reviews Physics3\(6\),pp\. 422–440\.External Links:ISSN 2522\-5820,[Document](https://dx.doi.org/10.1038/s42254-021-00314-5),[Link](https://doi.org/10.1038/s42254-021-00314-5)Cited by:[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p1.1)\.
- \[19\]A\. Gholami, K\. Keutzer, and G\. Biros\(2019\)ANODE: unconditionally accurate memory\-efficient gradients for neural ODEs\.Note:arXivExternal Links:1902\.10298,[Link](https://arxiv.org/abs/1902.10298)Cited by:[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p2.1)\.
- \[20\]J\. Zhuang, N\. Dvornek, X\. Li, S\. Tatikonda, X\. Papademetris, and J\. Duncan\(2020\)Adaptive checkpoint adjoint method for gradient estimation in neural ODE\.Note:arXivExternal Links:2006\.02493,[Link](https://arxiv.org/abs/2006.02493)Cited by:[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p2.1)\.
- \[21\]D\. Onken and L\. Ruthotto\(2020\)Discretize\-optimize vs\. optimize\-discretize for time\-series regression and continuous normalizing flows\.Note:arXivExternal Links:2005\.13420,[Link](https://arxiv.org/abs/2005.13420)Cited by:[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p2.1)\.
- \[22\]E\. M\. Turan and J\. Jäschke\(2022\)Multiple shooting for training neural differential equations on time series\.IEEE Control Systems Letters6,pp\. 1897–1902\.External Links:[Document](https://dx.doi.org/10.1109/LCSYS.2021.3135835)Cited by:[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p2.1)\.
- \[23\]E\. Roesch, C\. Rackauckas, and M\. P\. H\. Stumpf\(2021\)Collocation based training of neural ordinary differential equations\.Statistical Applications in Genetics and Molecular Biology20\(2\),pp\. 37–49\.External Links:[Document](https://dx.doi.org/10.1515/sagmb-2020-0025),[Link](https://doi.org/10.1515/sagmb-2020-0025)Cited by:[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p2.1)\.
- \[24\]J\. C\. Butcher\(2008\)Numerical methods for ordinary differential equations\.2 edition,John Wiley & Sons, Ltd\.External Links:ISBN 978\-0\-470\-72335\-7,[Document](https://dx.doi.org/10.1002/9780470753767)Cited by:[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p2.1)\.
- \[25\]E\. Hairer, C\. Lubich, and G\. Wanner\(2006\)Geometric numerical integration: structure\-preserving algorithms for ordinary differential equations\.2 edition,Springer Series in Computational Mathematics, Vol\.31,Springer Berlin, Heidelberg\.External Links:ISBN 978\-3\-540\-30663\-4,[Document](https://dx.doi.org/10.1007/3-540-30666-8)Cited by:[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p2.1)\.
- \[26\]T\. Shinbrot, C\. Grebogi, J\. Wisdom, and J\. A\. Yorke\(1992\)Chaos in a double pendulum\.American Journal of Physics60\(6\),pp\. 491–499\.External Links:[Document](https://dx.doi.org/10.1119/1.16860)Cited by:[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p2.1)\.
- \[27\]S\. H\. Strogatz\(2018\)Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering\.2 edition,CRC Press\.External Links:ISBN 978\-0\-8133\-4910\-7,[Link](https://www.biodyn.ro/course/literatura/Nonlinear_Dynamics_and_Chaos_2018_Steven_H._Strogatz.pdf)Cited by:[§2\.1](https://arxiv.org/html/2609.26972#S2.SS1.p2.1)\.
- \[28\]P\. Warden and D\. Situnayake\(2020\)TinyML: machine learning with TensorFlow Lite on Arduino and ultra\-low\-power microcontrollers\.1 edition,O’Reilly Media\.External Links:ISBN 978\-1\-4920\-5204\-3,[Link](https://www.oreilly.com/library/view/tinyml/9781492052036/)Cited by:[§2\.2](https://arxiv.org/html/2609.26972#S2.SS2.p1.1)\.
- \[29\]H\. Ren, D\. Anicic, and T\. A\. Runkler\(2021\)TinyOL: TinyML with online\-learning on microcontrollers\.In2021 International Joint Conference on Neural Networks \(IJCNN\),pp\. 1–8\.External Links:ISBN 978\-1\-6654\-3900\-8,[Document](https://dx.doi.org/10.1109/IJCNN52387.2021.9533927),[Link](https://doi.org/10.1109/IJCNN52387.2021.9533927)Cited by:[§2\.2](https://arxiv.org/html/2609.26972#S2.SS2.p1.1)\.
- \[30\]C\. Banbury, V\. J\. Reddi, P\. Torelli, J\. Holleman, N\. Jeffries, C\. Kiraly, P\. Montino, D\. Kanter, S\. Ahmed, D\. Pau, U\. Thakker, A\. Torrini, P\. Warden, J\. Cordaro, G\. Di Guglielmo, J\. Duarte, S\. Gibellini, V\. Parekh, H\. Tran, N\. Tran, W\. Niu, and X\. Xu\(2021\)MLPerf tiny benchmark\.Note:arXivExternal Links:[Document](https://dx.doi.org/10.48550/arXiv.2106.07597),2106\.07597,[Link](https://arxiv.org/abs/2106.07597)Cited by:[§2\.2](https://arxiv.org/html/2609.26972#S2.SS2.p1.1)\.
- \[31\]S\. Dhar, J\. Guo, J\. Liu, S\. Tripathi, U\. Kurup, and M\. Shah\(2021\)A survey of on\-device machine learning\.ACM Transactions on Internet of Things2\(3\)\.External Links:[Document](https://dx.doi.org/10.1145/3450494)Cited by:[§2\.2](https://arxiv.org/html/2609.26972#S2.SS2.p1.1)\.
- \[32\]H\. Cai, C\. Gan, L\. Zhu, and S\. Han\(2020\)TinyTL: reduce memory, not parameters for efficient on\-device learning\.InAdvances in Neural Information Processing Systems,Vol\.33,pp\. 11285–11297\.External Links:[Link](https://proceedings.neurips.cc/paper/2020/hash/81f7acabd411274fcf65ce2070ed568a-Abstract.html)Cited by:[§2\.2](https://arxiv.org/html/2609.26972#S2.SS2.p1.1)\.
- \[33\]J\. C\. Spall\(1992\)Multivariate stochastic approximation using a simultaneous perturbation gradient approximation\.IEEE Transactions on Automatic Control37\(3\),pp\. 332–341\.External Links:[Document](https://dx.doi.org/10.1109/9.119632)Cited by:[§2\.2](https://arxiv.org/html/2609.26972#S2.SS2.p1.1)\.
- \[34\]J\. C\. Spall\(1998\)An overview of the simultaneous perturbation method for efficient optimization\.Johns Hopkins APL Technical Digest19\(4\),pp\. 482–492\.External Links:[Link](https://www.jhuapl.edu/Content/techdigest/pdf/V19-N04/19-04-Spall.pdf)Cited by:[§2\.2](https://arxiv.org/html/2609.26972#S2.SS2.p1.1)\.
- \[35\]V\. Balakrishnan\(2002\)System identification: theory for the user \(second edition\): Lennart Ljung; Prentice\-Hall, Englewood Cliffs, NJ, 1999, ISBN 0\-13\-656695\-2\.Automatica38\(2\),pp\. 375–378\.External Links:ISSN 0005\-1098,[Document](https://dx.doi.org/10.1016/S0005-1098%2801%2900214-X),[Link](https://www.sciencedirect.com/science/article/pii/S000510980100214X)Cited by:[§2\.3](https://arxiv.org/html/2609.26972#S2.SS3.p1.1)\.
- \[36\]O\. Nelles\(2001\)Nonlinear system identification: from classical approaches to neural networks and fuzzy models\.1 edition,Springer Berlin, Heidelberg\.External Links:ISBN 978\-3\-540\-67369\-9,[Document](https://dx.doi.org/10.1007/978-3-662-04323-3)Cited by:[§2\.3](https://arxiv.org/html/2609.26972#S2.SS3.p1.1)\.
- \[37\]J\. M\. Varah\(1982\)A spline least squares method for numerical parameter estimation in differential equations\.SIAM Journal on Scientific and Statistical Computing3\(1\),pp\. 28–46\.External Links:[Document](https://dx.doi.org/10.1137/0903003)Cited by:[§2\.3](https://arxiv.org/html/2609.26972#S2.SS3.p1.1)\.
- \[38\]J\. O\. Ramsay, G\. Hooker, D\. Campbell, and J\. Cao\(2007\)Parameter estimation for differential equations: a generalized smoothing approach\.Journal of the Royal Statistical Society Series B: Statistical Methodology69\(5\),pp\. 741–796\.External Links:ISSN 1369\-7412,[Document](https://dx.doi.org/10.1111/j.1467-9868.2007.00610.x),[Link](https://doi.org/10.1111/j.1467-9868.2007.00610.x)Cited by:[§2\.3](https://arxiv.org/html/2609.26972#S2.SS3.p1.1)\.
- \[39\]M\. Schmidt and H\. Lipson\(2009\)Distilling free\-form natural laws from experimental data\.Science324\(5923\),pp\. 81–85\.External Links:[Document](https://dx.doi.org/10.1126/science.1165893),[Link](https://www.science.org/doi/abs/10.1126/science.1165893)Cited by:[§2\.3](https://arxiv.org/html/2609.26972#S2.SS3.p1.1)\.
- \[40\]S\. L\. Brunton, J\. L\. Proctor, and J\. N\. Kutz\(2016\)Discovering governing equations from data by sparse identification of nonlinear dynamical systems\.Proceedings of the National Academy of Sciences113\(15\),pp\. 3932–3937\.External Links:[Document](https://dx.doi.org/10.1073/pnas.1517384113),[Link](https://www.pnas.org/doi/abs/10.1073/pnas.1517384113)Cited by:[§2\.3](https://arxiv.org/html/2609.26972#S2.SS3.p1.1),[Table 1](https://arxiv.org/html/2609.26972#S2.T1.2.4.1)\.
- \[41\]S\. L\. Brunton and J\. N\. Kutz\(2022\)Data\-driven science and engineering: machine learning, dynamical systems, and control\.2 edition,Cambridge University Press,Cambridge\.External Links:ISBN 978\-1\-009\-09848\-9,[Document](https://dx.doi.org/10.1017/9781009089517),[Link](https://www.cambridge.org/highereducation/isbn/9781009089517)Cited by:[§2\.3](https://arxiv.org/html/2609.26972#S2.SS3.p1.1)\.
- \[42\]K\. Champion, B\. Lusch, J\. N\. Kutz, and S\. L\. Brunton\(2019\)Data\-driven discovery of coordinates and governing equations\.Proceedings of the National Academy of Sciences116\(45\),pp\. 22445–22451\.External Links:[Document](https://dx.doi.org/10.1073/pnas.1906995116),[Link](https://www.pnas.org/doi/abs/10.1073/pnas.1906995116)Cited by:[§2\.3](https://arxiv.org/html/2609.26972#S2.SS3.p1.1)\.
- \[43\]K\. Kaheman, J\. N\. Kutz, and S\. L\. Brunton\(2020\)SINDy\-PI: a robust algorithm for parallel implicit sparse identification of nonlinear dynamics\.Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences476\(2242\),pp\. 20200279\.External Links:ISSN 1471\-2946,[Document](https://dx.doi.org/10.1098/rspa.2020.0279),[Link](https://doi.org/10.1098/rspa.2020.0279)Cited by:[§2\.3](https://arxiv.org/html/2609.26972#S2.SS3.p1.1)\.
- \[44\]R\. Chartrand\(2011\)Numerical differentiation of noisy, nonsmooth data\.ISRN Applied Mathematics2011,pp\. 164564\.External Links:[Document](https://dx.doi.org/10.5402/2011/164564),[Link](https://onlinelibrary.wiley.com/doi/abs/10.5402/2011/164564)Cited by:[§2\.3](https://arxiv.org/html/2609.26972#S2.SS3.p1.1)\.
- \[45\]F\. van Breugel, J\. N\. Kutz, and B\. W\. Brunton\(2020\)Numerical differentiation of noisy data: a unifying multi\-objective optimization framework\.IEEE Access8,pp\. 196865–196877\.External Links:[Document](https://dx.doi.org/10.1109/ACCESS.2020.3034077)Cited by:[§2\.3](https://arxiv.org/html/2609.26972#S2.SS3.p1.1)\.
- \[46\]R\. W\. Schafer\(2011\)What is a Savitzky\-Golay filter? \[lecture notes\]\.IEEE Signal Processing Magazine28\(4\),pp\. 111–117\.External Links:[Document](https://dx.doi.org/10.1109/MSP.2011.941097)Cited by:[§2\.3](https://arxiv.org/html/2609.26972#S2.SS3.p1.1)\.
- \[47\]P\. J\. Huber\(1964\)Robust Estimation of a Location Parameter\.The Annals of Mathematical Statistics35\(1\),pp\. 73–101\.External Links:[Document](https://dx.doi.org/10.1214/aoms/1177703732),[Link](https://doi.org/10.1214/aoms/1177703732)Cited by:[item 3](https://arxiv.org/html/2609.26972#S2.I1.i3.p1.1)\.
- \[48\]A\. Griewank and A\. Walther\(2008\)Evaluating derivatives: principles and techniques of algorithmic differentiation\.2 edition,Other Titles in Applied Mathematics,Society for Industrial and Applied Mathematics\.External Links:ISBN 978\-0\-89871\-659\-7,[Document](https://dx.doi.org/10.1137/1.9780898717761)Cited by:[§2\.4](https://arxiv.org/html/2609.26972#S2.SS4.p1.1)\.
- \[49\]A\. G\. Baydin, B\. A\. Pearlmutter, A\. A\. Radul, and J\. M\. Siskind\(2018\)Automatic differentiation in machine learning: a survey\.Note:arXivExternal Links:1502\.05767v4,[Link](https://arxiv.org/abs/1502.05767v4)Cited by:[§2\.4](https://arxiv.org/html/2609.26972#S2.SS4.p1.1),[§3\.2](https://arxiv.org/html/2609.26972#S3.SS2.p3.4)\.
- \[50\]B\. A\. Pearlmutter\(1994\)Fast exact multiplication by the Hessian\.Neural Computation6\(1\),pp\. 147–160\.External Links:[Document](https://dx.doi.org/10.1162/neco.1994.6.1.147)Cited by:[§2\.4](https://arxiv.org/html/2609.26972#S2.SS4.p1.1),[§3\.2](https://arxiv.org/html/2609.26972#S3.SS2.p3.4)\.
- \[51\]W\. Squire and G\. Trapp\(1998\)Using complex variables to estimate derivatives of real functions\.SIAM Review40\(1\),pp\. 110–112\.External Links:[Document](https://dx.doi.org/10.1137/S003614459631241X)Cited by:[§2\.4](https://arxiv.org/html/2609.26972#S2.SS4.p1.1)\.
- \[52\]J\. R\. R\. A\. Martins, P\. Sturdza, and J\. J\. Alonso\(2003\)The complex\-step derivative approximation\.ACM Transactions on Mathematical Software29\(3\),pp\. 245–262\.External Links:ISSN 0098\-3500,[Document](https://dx.doi.org/10.1145/838250.838251),[Link](https://doi.org/10.1145/838250.838251)Cited by:[§2\.4](https://arxiv.org/html/2609.26972#S2.SS4.p1.1)\.
- \[53\]N\. K\. Sinha and M\. P\. Griscik\(1971\)A stochastic approximation method\.IEEE Transactions on Systems, Man, and CyberneticsSMC\-1\(4\),pp\. 338–344\.External Links:[Document](https://dx.doi.org/10.1109/TSMC.1971.4308316)Cited by:[§2\.4](https://arxiv.org/html/2609.26972#S2.SS4.p1.1)\.
- \[54\]I\. Loshchilov and F\. Hutter\(2017\)SGDR: stochastic gradient descent with warm restarts\.Note:arXivExternal Links:1608\.03983v5,[Link](https://arxiv.org/abs/1608.03983v5)Cited by:[§2\.4](https://arxiv.org/html/2609.26972#S2.SS4.p1.1)\.
- \[55\]P\. Izmailov, D\. Podoprikhin, T\. Garipov, D\. Vetrov, and A\. G\. Wilson\(2019\)Averaging weights leads to wider optima and better generalization\.Note:arXivExternal Links:1803\.05407v3,[Link](https://arxiv.org/abs/1803.05407v3)Cited by:[§2\.4](https://arxiv.org/html/2609.26972#S2.SS4.p1.1)\.
- \[56\]J\. S\. Vitter\(1985\)Random sampling with a reservoir\.ACM Transactions on Mathematical Software11\(1\),pp\. 37–57\.External Links:ISSN 0098\-3500,[Document](https://dx.doi.org/10.1145/3147.3165)Cited by:[§2\.4](https://arxiv.org/html/2609.26972#S2.SS4.p1.1)\.

## Appendix ADerivations

### A\.1Reverse\-over\-JVP Derivation for Second\-Order Jet Gradients

The second\-order derivative loss term is given by:

ℒ2=λ2​∥r2∥22,\\mathcal\{L\}\_\{2\}=\\lambda\_\{2\}\\lVert r\_\{2\}\\rVert\_\{2\}^\{2\},\(28\)where residualr2r\_\{2\}and predicted second derivativeGθG\_\{\\theta\}are defined as:

r2=x¨^−Gθ​\(x^\)s2,r\_\{2\}=\\frac\{\\hat\{\\ddot\{x\}\}\-G\_\{\\theta\}\(\\hat\{x\}\)\}\{s\_\{2\}\},\(29\)Gθ​\(x^\)=JFθ​\(x^\)​Fθ​\(x^\)\.G\_\{\\theta\}\(\\hat\{x\}\)=J\_\{F\_\{\\theta\}\}\(\\hat\{x\}\)\\,F\_\{\\theta\}\(\\hat\{x\}\)\.\(30\)Differentiatingℒ2\\mathcal\{L\}\_\{2\}with respect to parametersθ\\thetayields:

d​ℒ2=−2​λ2​\(r2s2\)⊤​d​Gθ,\\mathrm\{d\}\\mathcal\{L\}\_\{2\}=\-2\\lambda\_\{2\}\\left\(\\frac\{r\_\{2\}\}\{s\_\{2\}\}\\right\)^\{\\\!\\top\}\\mathrm\{d\}G\_\{\\theta\},\(31\)where the differentiald​Gθ\\mathrm\{d\}G\_\{\\theta\}expands to:

d​Gθ=\(d​JNθ\)​Fθ​\(x^\)\+JFθ​\(x^\)​d​Nθ\.\\mathrm\{d\}G\_\{\\theta\}=\(\\mathrm\{d\}J\_\{N\_\{\\theta\}\}\)\\,F\_\{\\theta\}\(\\hat\{x\}\)\+J\_\{F\_\{\\theta\}\}\(\\hat\{x\}\)\\,\\mathrm\{d\}N\_\{\\theta\}\.\(32\)Substituting weight vectorw=−2​\(r2/s2\)w=\-2\(r\_\{2\}/s\_\{2\}\), we separate the gradient into two parts:

d​ℒ2=λ2​\[w⊤​\(d​JNθ\)​v\+w⊤​JFθ​\(x^\)​d​Nθ\],\\mathrm\{d\}\\mathcal\{L\}\_\{2\}=\\lambda\_\{2\}\\left\[w^\{\\\!\\top\}\(\\mathrm\{d\}J\_\{N\_\{\\theta\}\}\)v\+w^\{\\\!\\top\}J\_\{F\_\{\\theta\}\}\(\\hat\{x\}\)\\,\\mathrm\{d\}N\_\{\\theta\}\\right\],\(33\)wherev=Fθ​\(x^\)v=F\_\{\\theta\}\(\\hat\{x\}\)\. Using the vectorization identitya⊤​X​b=vec​\(X\)⊤​\(b⊗a\)a^\{\\\!\\top\}Xb=\\mathrm\{vec\}\(X\)^\{\\\!\\top\}\(b\\otimes a\), the first term transforms to:

w⊤​\(d​JNθ\)​v=vec​\(d​JNθ\)⊤​\(v⊗w\)\.w^\{\\\!\\top\}\(\\mathrm\{d\}J\_\{N\_\{\\theta\}\}\)v=\\mathrm\{vec\}\(\\mathrm\{d\}J\_\{N\_\{\\theta\}\}\)^\{\\\!\\top\}\(v\\otimes w\)\.\(34\)This allows computing parameter gradients by reverse\-mode differentiation over the Jacobian–vector productJNθ​vJ\_\{N\_\{\\theta\}\}v\(the forward\-mode tangent pass\) without explicitly forming higher\-order Jacobian matrices\.

### A\.2Noise Estimator Variance Identity

Let a local signal window be represented as𝐲=𝐬\+𝜼\\mathbf\{y\}=\\mathbf\{s\}\+\\boldsymbol\{\\eta\}, where noise𝜼∼𝒩⁡\(0,σ2​I\)\\boldsymbol\{\\eta\}\\sim\\mathcal\{N\}\(0,\\sigma^\{2\}I\)and𝐬=A​𝐜\\mathbf\{s\}=A\\mathbf\{c\}represents an exact polynomial sequence\. Because Savitzky–Golay smoothing tapsD0D\_\{0\}preserve polynomials at the center indexμ\\mu, we haveD0⊤​A​𝐜=sμD\_\{0\}^\{\\\!\\top\}A\\mathbf\{c\}=s\_\{\\mu\}\. The scalar fit residual is:

e=yμ−D0⊤​𝐲=ημ−D0⊤​𝜼\.e=y\_\{\\mu\}\-D\_\{0\}^\{\\\!\\top\}\\mathbf\{y\}=\\eta\_\{\\mu\}\-D\_\{0\}^\{\\\!\\top\}\\boldsymbol\{\\eta\}\.\(35\)Taking the variance yields:

Var⁡\(e\)=Var⁡\(ημ\)−2​Cov​\(ημ,D0⊤​𝜼\)\+Var⁡\(D0⊤​𝜼\)\.\\mathrm\{Var\}\(e\)=\\mathrm\{Var\}\(\\eta\_\{\\mu\}\)\-2\\mathrm\{Cov\}\(\\eta\_\{\\mu\},D\_\{0\}^\{\\\!\\top\}\\boldsymbol\{\\eta\}\)\+\\mathrm\{Var\}\(D\_\{0\}^\{\\\!\\top\}\\boldsymbol\{\\eta\}\)\.\(36\)Substituting component variancesVar⁡\(ημ\)=σ2\\mathrm\{Var\}\(\\eta\_\{\\mu\}\)=\\sigma^\{2\},Cov⁡\(ημ,D0⊤​𝜼\)=σ2​\[D0\]μ\\mathrm\{Cov\}\(\\eta\_\{\\mu\},D\_\{0\}^\{\\\!\\top\}\\boldsymbol\{\\eta\}\)=\\sigma^\{2\}\[D\_\{0\}\]\_\{\\mu\}, andVar⁡\(D0⊤​𝜼\)=σ2​∥D0∥22\\mathrm\{Var\}\(D\_\{0\}^\{\\\!\\top\}\\boldsymbol\{\\eta\}\)=\\sigma^\{2\}\\lVert D\_\{0\}\\rVert\_\{2\}^\{2\}gives:

Var⁡\(e\)=σ2​\(1−2​\[D0\]μ\+∥D0∥22\)\.\\mathrm\{Var\}\(e\)=\\sigma^\{2\}\\left\(1\-2\[D\_\{0\}\]\_\{\\mu\}\+\\lVert D\_\{0\}\\rVert\_\{2\}^\{2\}\\right\)\.\(37\)Inverting this expression provides the noise variance estimator in \([23](https://arxiv.org/html/2609.26972#S3.E23)\)\.

### A\.3Derivation of Damped Single Pendulum Dynamics

Consider a pendulum of massmmattached to a rigid rod of lengthLLpivoted at the origin\(0,0\)\(0,0\), with angular displacementθ\\thetameasured relative to the downward vertical axis\. The Cartesian position vector of the pendulum bob isr→=\(L​sin⁡θ\)​i^\+\(−L​cos⁡θ\)​j^\\vec\{r\}=\(L\\sin\\theta\)\\hat\{i\}\+\(\-L\\cos\\theta\)\\hat\{j\}\. Taking the time derivative gives the velocity vectorv→=x˙​i^\+y˙​j^\\vec\{v\}=\\dot\{x\}\\hat\{i\}\+\\dot\{y\}\\hat\{j\}, where application of the chain rule yields:

x˙\\displaystyle\\dot\{x\}=dd​t​\(L​sin⁡θ\)=L​θ˙​cos⁡θ,\\displaystyle=\\frac\{d\}\{dt\}\\left\(L\\sin\\theta\\right\)=L\\dot\{\\theta\}\\cos\\theta,\(38\)y˙\\displaystyle\\dot\{y\}=dd​t​\(−L​cos⁡θ\)=L​θ˙​sin⁡θ\.\\displaystyle=\\frac\{d\}\{dt\}\\left\(\-L\\cos\\theta\\right\)=L\\dot\{\\theta\}\\sin\\theta\.\(39\)
The magnitude squared of the velocity vector is‖v→‖2=x˙2\+y˙2\\\|\\vec\{v\}\\\|^\{2\}=\\dot\{x\}^\{2\}\+\\dot\{y\}^\{2\}\. Substituting these velocity components into the kinetic energy formulationTTand applying the identitycos2⁡θ\+sin2⁡θ=1\\cos^\{2\}\\theta\+\\sin^\{2\}\\theta=1gives:

T=12​m​\(x˙2\+y˙2\)=12​m​\(L2​θ˙2​cos2⁡θ\+L2​θ˙2​sin2⁡θ\)=12​m​L2​θ˙2\.T=\\frac\{1\}\{2\}m\(\\dot\{x\}^\{2\}\+\\dot\{y\}^\{2\}\)=\\frac\{1\}\{2\}m\\left\(L^\{2\}\\dot\{\\theta\}^\{2\}\\cos^\{2\}\\theta\+L^\{2\}\\dot\{\\theta\}^\{2\}\\sin^\{2\}\\theta\\right\)=\\frac\{1\}\{2\}mL^\{2\}\\dot\{\\theta\}^\{2\}\.\(40\)Takingy=0y=0as the reference plane, the potential energy isV=m​g​y=−m​g​L​cos⁡θV=mgy=\-mgL\\cos\\theta\. The system LagrangianLg=T−VL\_\{g\}=T\-Vis thus formulated as:

Lg=12​m​L2​θ˙2\+m​g​L​cos⁡θ\.L\_\{g\}=\\frac\{1\}\{2\}mL^\{2\}\\dot\{\\theta\}^\{2\}\+mgL\\cos\\theta\.\(41\)
To model mechanical dissipation, we introduce a generalized non\-conservative viscous damping forceQnc=−b​θ˙Q\_\{\\text\{nc\}\}=\-b\\dot\{\\theta\}, wherebbis the rotational damping coefficient\. Applying the generalized Euler\-Lagrange equation with dissipation:

dd​t​\(∂Lg∂θ˙\)−∂Lg∂θ=Qnc\.\\frac\{d\}\{dt\}\\left\(\\frac\{\\partial L\_\{g\}\}\{\\partial\\dot\{\\theta\}\}\\right\)\-\\frac\{\\partial L\_\{g\}\}\{\\partial\\theta\}=Q\_\{\\text\{nc\}\}\.\(42\)
Evaluating the partial derivatives yields∂Lg∂θ˙=m​L2​θ˙⟹dd​t​\(∂Lg∂θ˙\)=m​L2​θ¨\\frac\{\\partial L\_\{g\}\}\{\\partial\\dot\{\\theta\}\}=mL^\{2\}\\dot\{\\theta\}\\implies\\frac\{d\}\{dt\}\\left\(\\frac\{\\partial L\_\{g\}\}\{\\partial\\dot\{\\theta\}\}\\right\)=mL^\{2\}\\ddot\{\\theta\}and∂Lg∂θ=−m​g​L​sin⁡θ\\frac\{\\partial L\_\{g\}\}\{\\partial\\theta\}=\-mgL\\sin\\theta\. Substituting these terms into the Euler\-Lagrange equation produces:

m​L2​θ¨−\(−m​g​L​sin⁡θ\)=−b​θ˙⟹θ¨\+gL​sin⁡θ=−bm​L2​θ˙\.mL^\{2\}\\ddot\{\\theta\}\-\(\-mgL\\sin\\theta\)=\-b\\dot\{\\theta\}\\implies\\ddot\{\\theta\}\+\\frac\{g\}\{L\}\\sin\\theta=\-\\frac\{b\}\{mL^\{2\}\}\\dot\{\\theta\}\.\(43\)
Defining state variablesx1=θx\_\{1\}=\\thetaandx2=ω=θ˙x\_\{2\}=\\omega=\\dot\{\\theta\}, alongside the normalized damping factorc=bm​L2c=\\frac\{b\}\{mL^\{2\}\}, reduces the second\-order differential equation to the first\-order autonomous state\-space system:

θ˙\\displaystyle\\dot\{\\theta\}=ω,\\displaystyle=\\omega,\(44\)ω˙\\displaystyle\\dot\{\\omega\}=−gL​sin⁡θ−c​ω\.\\displaystyle=\-\\frac\{g\}\{L\}\\sin\\theta\-c\\omega\.\(45\)

### A\.4Derivation of Damped Double Pendulum Dynamics

Consider a double pendulum composed of two point massesm1m\_\{1\}andm2m\_\{2\}suspended by rigid, massless rods of lengthsl1l\_\{1\}andl2l\_\{2\}, respectively\. Letθ1\\theta\_\{1\}andθ2\\theta\_\{2\}denote the angular displacements of the upper and lower links relative to the downward vertical axis\.

The Cartesian coordinates ofm1m\_\{1\}andm2m\_\{2\}are given by:

x1\\displaystyle x\_\{1\}=l1sinθ1,\\displaystyle=l\_\{1\}\\sin\\theta\_\{1\},y1\\displaystyle y\_\{1\}=−l1cosθ1\\displaystyle=\-l\_\{1\}\\cos\\theta\_\{1\}\(46\)x2\\displaystyle x\_\{2\}=l1sinθ1\+l2sinθ2,\\displaystyle=l\_\{1\}\\sin\\theta\_\{1\}\+l\_\{2\}\\sin\\theta\_\{2\},y2\\displaystyle y\_\{2\}=−l1cosθ1−l2cosθ2\\displaystyle=\-l\_\{1\}\\cos\\theta\_\{1\}\-l\_\{2\}\\cos\\theta\_\{2\}\(47\)
Differentiating with respect to timettyields the linear velocity components\(x˙1,y˙1\)\(\\dot\{x\}\_\{1\},\\dot\{y\}\_\{1\}\)and\(x˙2,y˙2\)\(\\dot\{x\}\_\{2\},\\dot\{y\}\_\{2\}\)\. Substituting these components into the total kinetic energy expressionT=12​m1​\(x˙12\+y˙12\)\+12​m2​\(x˙22\+y˙22\)T=\\frac\{1\}\{2\}m\_\{1\}\(\\dot\{x\}\_\{1\}^\{2\}\+\\dot\{y\}\_\{1\}^\{2\}\)\+\\frac\{1\}\{2\}m\_\{2\}\(\\dot\{x\}\_\{2\}^\{2\}\+\\dot\{y\}\_\{2\}^\{2\}\)and applying trigonometric angle addition identities produces:

T=12​\(m1\+m2\)​l12​θ˙12\+12​m2​l22​θ˙22\+m2​l1​l2​θ˙1​θ˙2​cos⁡\(θ1−θ2\)T=\\frac\{1\}\{2\}\(m\_\{1\}\+m\_\{2\}\)l\_\{1\}^\{2\}\\dot\{\\theta\}\_\{1\}^\{2\}\+\\frac\{1\}\{2\}m\_\{2\}l\_\{2\}^\{2\}\\dot\{\\theta\}\_\{2\}^\{2\}\+m\_\{2\}l\_\{1\}l\_\{2\}\\dot\{\\theta\}\_\{1\}\\dot\{\\theta\}\_\{2\}\\cos\(\\theta\_\{1\}\-\\theta\_\{2\}\)\(48\)
Taking the fixed pivot at the origin\(0,0\)\(0,0\)as the potential energy reference plane \(y=0y=0\), the potential energyVVis:

V=m1gy1\+m2gy2=−\(m1\+m2\)gl1cosθ1−m2gl2cosθ2V=m\_\{1\}gy\_\{1\}\+m\_\{2\}gy\_\{2\}=\-\(m\_\{1\}\+m\_\{2\}\)gl\_\{1\}\\cos\\theta\_\{1\}\-m\_\{2\}gl\_\{2\}\\cos\\theta\_\{2\}\(49\)
Formulating the system LagrangianLg=T−VL\_\{g\}=T\-Vand introducing Rayleigh dissipation forcesQnc,1=−c1​θ˙1Q\_\{\\text\{nc\},1\}=\-c\_\{1\}\\dot\{\\theta\}\_\{1\}andQnc,2=−c2​θ˙2Q\_\{\\text\{nc\},2\}=\-c\_\{2\}\\dot\{\\theta\}\_\{2\}, we evaluate the coupled Euler\-Lagrange equationsdd​t​\(∂Lg∂θ˙i\)−∂Lg∂θi=Qnc,i\\frac\{d\}\{dt\}\\left\(\\frac\{\\partial L\_\{g\}\}\{\\partial\\dot\{\\theta\}\_\{i\}\}\\right\)\-\\frac\{\\partial L\_\{g\}\}\{\\partial\\theta\_\{i\}\}=Q\_\{\\text\{nc\},i\}fori∈\{1,2\}i\\in\\\{1,2\\\}\.

Defining angular velocitiesω1=θ˙1\\omega\_\{1\}=\\dot\{\\theta\}\_\{1\},ω2=θ˙2\\omega\_\{2\}=\\dot\{\\theta\}\_\{2\}, and the angular differenceδ=θ1−θ2\\delta=\\theta\_\{1\}\-\\theta\_\{2\}, the resulting system can be expressed as a linear matrix equation in angular accelerations\(ω˙1,ω˙2\)\(\\dot\{\\omega\}\_\{1\},\\dot\{\\omega\}\_\{2\}\):

\[\(m1\+m2\)​l1m2​l2​cos⁡δm2​l1​cos⁡δm2​l2\]​\[ω˙1ω˙2\]=\[−m2l2ω22sinδ−\(m1\+m2\)gsinθ1−c1ω1m2l1ω12sinδ−m2gsinθ2−c2ω2\]\\begin\{bmatrix\}\(m\_\{1\}\+m\_\{2\}\)l\_\{1\}&m\_\{2\}l\_\{2\}\\cos\\delta\\\\ m\_\{2\}l\_\{1\}\\cos\\delta&m\_\{2\}l\_\{2\}\\end\{bmatrix\}\\begin\{bmatrix\}\\dot\{\\omega\}\_\{1\}\\\\ \\dot\{\\omega\}\_\{2\}\\end\{bmatrix\}=\\begin\{bmatrix\}\-m\_\{2\}l\_\{2\}\\omega\_\{2\}^\{2\}\\sin\\delta\-\(m\_\{1\}\+m\_\{2\}\)g\\sin\\theta\_\{1\}\-c\_\{1\}\\omega\_\{1\}\\\\ m\_\{2\}l\_\{1\}\\omega\_\{1\}^\{2\}\\sin\\delta\-m\_\{2\}g\\sin\\theta\_\{2\}\-c\_\{2\}\\omega\_\{2\}\\end\{bmatrix\}\(50\)
Solving the2×22\\times 2linear system via Cramer’s Rule requires the determinant of the mass matrixMM:

det\(M\)=m2​l1​l2​\(m1\+m2−m2​cos2⁡δ\)\\det\(M\)=m\_\{2\}l\_\{1\}l\_\{2\}\\left\(m\_\{1\}\+m\_\{2\}\-m\_\{2\}\\cos^\{2\}\\delta\\right\)\(51\)
Using the double\-angle identitycos2⁡δ=1\+cos⁡\(2​δ\)2\\cos^\{2\}\\delta=\\frac\{1\+\\cos\(2\\delta\)\}\{2\}, we define the common denominator parameterDD:

D=2​\(m1\+m2−m2​cos2⁡δ\)=2​m1\+m2−m2​cos⁡\(2​δ\)D=2\(m\_\{1\}\+m\_\{2\}\-m\_\{2\}\\cos^\{2\}\\delta\)=2m\_\{1\}\+m\_\{2\}\-m\_\{2\}\\cos\(2\\delta\)\(52\)
Inverting the matrix system and isolatingω˙1\\dot\{\\omega\}\_\{1\}andω˙2\\dot\{\\omega\}\_\{2\}yields the explicit state\-space system equations:

θ˙1\\displaystyle\\dot\{\\theta\}\_\{1\}=ω1\\displaystyle=\\omega\_\{1\}\(53\)θ˙2\\displaystyle\\dot\{\\theta\}\_\{2\}=ω2\\displaystyle=\\omega\_\{2\}\(54\)ω˙1\\displaystyle\\dot\{\\omega\}\_\{1\}=−g\(2m1\+m2\)sinθ1−m2gsin\(θ1−2θ2\)−2sinδm2\(ω22l2\+ω12l1cosδ\)l1​D−c1​ω1\\displaystyle=\\frac\{\-g\(2m\_\{1\}\+m\_\{2\}\)\\sin\\theta\_\{1\}\-m\_\{2\}g\\sin\(\\theta\_\{1\}\-2\\theta\_\{2\}\)\-2\\sin\\delta\\,m\_\{2\}\\left\(\\omega\_\{2\}^\{2\}l\_\{2\}\+\\omega\_\{1\}^\{2\}l\_\{1\}\\cos\\delta\\right\)\}\{l\_\{1\}D\}\-c\_\{1\}\\omega\_\{1\}\(55\)ω˙2\\displaystyle\\dot\{\\omega\}\_\{2\}=2sinδ\(ω12l1\(m1\+m2\)\+g\(m1\+m2\)cosθ1\+ω22l2m2cosδ\)l2​D−c2​ω2\\displaystyle=\\frac\{2\\sin\\delta\\left\(\\omega\_\{1\}^\{2\}l\_\{1\}\(m\_\{1\}\+m\_\{2\}\)\+g\(m\_\{1\}\+m\_\{2\}\)\\cos\\theta\_\{1\}\+\\omega\_\{2\}^\{2\}l\_\{2\}m\_\{2\}\\cos\\delta\\right\)\}\{l\_\{2\}D\}\-c\_\{2\}\\omega\_\{2\}\(56\)

## Appendix BHyperparameter Reference

Table[5](https://arxiv.org/html/2609.26972#A2.T5)lists all configuration parameters used across single and double pendulum experiments\.

Table 5:System parameters and configuration settings\.
## Appendix CAnalysis of High\-Velocity Error Outliers

Experimental runs on the double pendulum identified localized error increases during rapid state transitions \(m=24m=24\)\. In failing seeds, training loss remained low \(5\.9×10−35\.9\{\\times\}10^\{\-3\}\) while field errors exceeded0\.110\.11\. Evaluating error distributions across state space showed that localized discrepancies occurred primarily at high angular velocities \(e\.g\.,x=\[−0\.73,−0\.30,4\.67,−7\.80\]x=\[\-0\.73,\-0\.30,4\.67,\-7\.80\]\), where high derivative magnitudes caused polynomial filter approximation errors\.

Because neural networks can fit biased targets tightly, standard robust loss functions \(such as Huber losses\) do not automatically reject these points\. Quality gating via polynomial misfit verification \([27](https://arxiv.org/html/2609.26972#S3.E27)\) detects these violations directly from raw sensor windows, reducing peak field errors from0\.1350\.135to0\.0250\.025\.

## Appendix DExtended Per\-Seed Results

Table[6](https://arxiv.org/html/2609.26972#A4.T6)provides detailed numerical results across individual random seeds\.

Table 6:Individual per\-seed relative field error values across test configurations\.
## Appendix ESupplementary Gradient and Estimator Validation Plots

This section contains supplementary validation plots for numerical gradient checks, Savitzky–Golay filter characteristics, and baseline sensitivity evaluations\.

Figure 10:Validation of analytic jet matching gradients against complex\-step differentiation across 32 test conditions\.Figure 11:Validation of baseline discrete adjoint gradients against complex\-step differentiation across 14 test conditions\.Figure 12:Savitzky–Golay derivative estimation error across effective step sizesΔ​teff\\Delta t\_\{\\mathrm\{eff\}\}and noise levels\.Figure 13:Effect of observation sampling rate on baseline adjoint field error\.