Neural Ordinary Differential Equations

• Neural Ordinary Differential Equations are continuous-depth models where hidden states evolve as ODE solutions rather than through discrete layers.
• They integrate adaptive numerical solvers into training, enabling constant-memory backpropagation and explicit trade-offs between accuracy and computational cost.
• Applications span from continuous-depth residual networks and latent variable models to continuous normalizing flows and structured scientific simulations.

Neural ordinary differential equations (Neural ODEs) are continuous-depth neural models in which the hidden state evolves as the solution of an ordinary differential equation rather than through a finite, explicitly enumerated stack of layers. In the canonical formulation, the state $\mathbf{h}(t)\in\mathbb{R}^D$ satisfies

$$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$

with output computed by an ODE solver over a time interval, so that evaluation and training become inseparable from numerical integration. This construction replaces discrete depth by a continuous variable, allows black-box solver usage, yields constant-memory backpropagation in theory, and exposes explicit trade-offs among accuracy, runtime, and solver tolerance. Since its introduction, the framework has expanded from continuous-depth residual networks to latent-variable models, continuous normalizing flows, scientific machine learning, uncertainty-aware modeling, and structured dynamical-system learning (Chen et al., 2018, Ruthotto, 2024).

1. Continuous-depth formulation and dynamical interpretation

At the center of a Neural ODE is the time-dependent vector field $f$, typically represented by a neural network. Given an initial condition $\mathbf{h}(t_0)$, the terminal state is

$$\mathbf{h}(t_1)=\mathbf{h}(t_0)+\int_{t_0}^{t_1} f\bigl(\mathbf{h}(t),t;\theta\bigr)\,dt =\mathrm{ODESolve}\bigl(\mathbf{h}(t_0),f,t_0,t_1,\theta\bigr).$$

The black-box integrator adaptively places as many function evaluations as needed to meet user-specified error tolerances. In the continuous-depth reading of residual architectures, the ResNet update $x^{k+1}=x^k+f(x^k;\theta^k)$ is interpreted as a forward Euler discretization of an ODE, and the continuous-depth limit is recovered as the step size tends to zero (Chen et al., 2018, Ruthotto, 2024).

This dynamical perspective immediately imports classical ODE concepts. Under a Lipschitz condition on $f(\cdot,t;\theta)$, Picard–Lindelöf yields existence and uniqueness of solutions, and the learned network becomes the time-$1$ flow map $\phi_\theta:x\mapsto h(1)$. The survey literature emphasizes that continuity and differentiability of $\phi_\theta$ with respect to both state and parameters follow from standard dependence-on-parameters results, allowing analysis through optimal control, transport equations, and continuity equations rather than only through layerwise neural-network heuristics (Ruthotto, 2024).

Solver choice is part of the model definition rather than a purely external implementation detail. Explicit solvers such as Runge–Kutta methods are simple and often fast per step; implicit solvers such as backward Euler or BDF are preferable for stiff dynamics. The original Neural ODE experiments used LSODE/VODE, Dormand–Prince, and related methods, with problem-dependent tolerances such as $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$0 for supervised classification, $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$1 for continuous normalizing flows, and a default near $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$2 for latent ODEs (Chen et al., 2018).

2. Adjoint sensitivity, memory behavior, and the continuous–discrete interface

Training requires gradients through the ODE solution operator. For a scalar terminal loss $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$3, the adjoint state

$$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$4

satisfies the backward ODE

$$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$5

The parameter gradient is then recovered by integrating

$$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$6

This adjoint sensitivity method enables backpropagation through any ODE solver without access to the solver’s internal operations and, because intermediate RK/Euler states need not be stored, gives $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$7 memory cost in the number of implied layers (Chen et al., 2018, Ruthotto, 2024).

The computational trade-off is explicit. Lower absolute and relative tolerances generally increase the number of function evaluations and raise runtime, while higher tolerances reduce cost at the expense of numerical precision. Adaptive variable-step solvers therefore allocate more computation to harder regions of the trajectory and less to smoother ones. The original experiments further report that backward NFE is often $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$8 forward NFE, so adjoint backpropagation can be faster than naïve differentiation through all solver steps (Chen et al., 2018).

The relation to discrete residual networks is more delicate than the informal Euler analogy suggests. A precise analysis of ResNets as discretizations of Neural ODEs shows that the distance between the discrete hidden-state trajectory and the corresponding continuous solution need not vanish with depth if the residual functions are not smooth with depth. Under linear residual functions and small enough initial loss, gradient descent preserves the required smoothness and induces convergence toward a limit Neural ODE at rate $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$9, uniformly with depth and optimization time. The same analysis shows that a memory-free discrete adjoint succeeds at large depth when the residual functions are Lipschitz with the input, and that Heun’s method provides better gradient estimation than Euler when the residual functions are smooth with depth (Sander et al., 2022).

Practical caveats remain. Mini-batching is implemented by concatenating states, which may inflate solver cost per batch; integrating the adjoint backward may accumulate extra numerical error; and tolerance choice must be tuned because it directly shifts the computability–precision trade-off. Checkpointing can ameliorate reversibility errors, though it was reported as unnecessary in the original work (Chen et al., 2018).

3. Canonical model classes and representative applications

The first canonical use of the framework is the continuous-depth residual network. Residual blocks of the form $f$0 are replaced by integrating $f$1 over a depth interval. On MNIST, ODE-Nets match ResNet accuracy at about $f$2 while maintaining $f$3 memory in depth. Their number of function evaluations adapts per example and increases during training, reflecting a compute-for-capacity trade-off (Chen et al., 2018).

A second canonical model class is the continuous-time latent variable model. In the latent ODE construction, one samples $f$4, evolves the latent path by solving $f$5, and emits observations from $f$6. Approximate inference is performed with an ODE-driven variational autoencoder, typically using an RNN encoder for $f$7. The model handles arbitrarily spaced timestamps and can incorporate an event-time likelihood through an inhomogeneous Poisson intensity $f$8. In the reported irregularly sampled spiral extrapolation experiment, latent ODEs achieved RMSE about $f$9 versus about $\mathbf{h}(t_0)$0 for RNN baselines (Chen et al., 2018).

The third canonical class is the continuous normalizing flow (CNF). If $\mathbf{h}(t_0)$1, the log-density evolves via the instantaneous change-of-variables formula

$\mathbf{h}(t_0)$2

This replaces the log-determinant of a full Jacobian by a Jacobian trace, which can be estimated in $\mathbf{h}(t_0)$3 time via Hutchinson’s trace estimator. CNFs therefore support exact likelihood training by integrating the augmented state for $\mathbf{h}(t_0)$4 and $\mathbf{h}(t_0)$5, and they avoid partitioning or ordering the data dimensions. The original experiments report better density matching than stacked planar flows of the same size (Chen et al., 2018).

Subsequent application work generalized the framework beyond image classification and density modeling. In semantic segmentation, SegNode reinterprets the four-branch, multi-resolution residual HRNetV2 module as the vector field inside an ODEBlock, drops explicit time dependence, and ties parameters across continuous depth. On CamVid with input size $\mathbf{h}(t_0)$6 and batch size $\mathbf{h}(t_0)$7, the discrete baseline with $\mathbf{h}(t_0)$8 M parameters required $\mathbf{h}(t_0)$9 GB memory for training and $$\mathbf{h}(t_1)=\mathbf{h}(t_0)+\int_{t_0}^{t_1} f\bigl(\mathbf{h}(t),t;\theta\bigr)\,dt =\mathrm{ODESolve}\bigl(\mathbf{h}(t_0),f,t_0,t_1,\theta\bigr).$$0 GB for single-image inference, whereas SegNode with $$\mathbf{h}(t_1)=\mathbf{h}(t_0)+\int_{t_0}^{t_1} f\bigl(\mathbf{h}(t),t;\theta\bigr)\,dt =\mathrm{ODESolve}\bigl(\mathbf{h}(t_0),f,t_0,t_1,\theta\bigr).$$1 M parameters required $$\mathbf{h}(t_1)=\mathbf{h}(t_0)+\int_{t_0}^{t_1} f\bigl(\mathbf{h}(t),t;\theta\bigr)\,dt =\mathrm{ODESolve}\bigl(\mathbf{h}(t_0),f,t_0,t_1,\theta\bigr).$$2 GB for training and $$\mathbf{h}(t_1)=\mathbf{h}(t_0)+\int_{t_0}^{t_1} f\bigl(\mathbf{h}(t),t;\theta\bigr)\,dt =\mathrm{ODESolve}\bigl(\mathbf{h}(t_0),f,t_0,t_1,\theta\bigr).$$3 GB for inference. Reported benchmark performance includes $$\mathbf{h}(t_1)=\mathbf{h}(t_0)+\int_{t_0}^{t_1} f\bigl(\mathbf{h}(t),t;\theta\bigr)\,dt =\mathrm{ODESolve}\bigl(\mathbf{h}(t_0),f,t_0,t_1,\theta\bigr).$$4 mIoU on Cityscapes, $$\mathbf{h}(t_1)=\mathbf{h}(t_0)+\int_{t_0}^{t_1} f\bigl(\mathbf{h}(t),t;\theta\bigr)\,dt =\mathrm{ODESolve}\bigl(\mathbf{h}(t_0),f,t_0,t_1,\theta\bigr).$$5 on CamVid, $$\mathbf{h}(t_1)=\mathbf{h}(t_0)+\int_{t_0}^{t_1} f\bigl(\mathbf{h}(t),t;\theta\bigr)\,dt =\mathrm{ODESolve}\bigl(\mathbf{h}(t_0),f,t_0,t_1,\theta\bigr).$$6 on LIP, and $$\mathbf{h}(t_1)=\mathbf{h}(t_0)+\int_{t_0}^{t_1} f\bigl(\mathbf{h}(t),t;\theta\bigr)\,dt =\mathrm{ODESolve}\bigl(\mathbf{h}(t_0),f,t_0,t_1,\theta\bigr).$$7 on PASCAL-Context, with SegNode + OCR refinement reaching $$\mathbf{h}(t_1)=\mathbf{h}(t_0)+\int_{t_0}^{t_1} f\bigl(\mathbf{h}(t),t;\theta\bigr)\,dt =\mathrm{ODESolve}\bigl(\mathbf{h}(t_0),f,t_0,t_1,\theta\bigr).$$8, $$\mathbf{h}(t_1)=\mathbf{h}(t_0)+\int_{t_0}^{t_1} f\bigl(\mathbf{h}(t),t;\theta\bigr)\,dt =\mathrm{ODESolve}\bigl(\mathbf{h}(t_0),f,t_0,t_1,\theta\bigr).$$9, $x^{k+1}=x^k+f(x^k;\theta^k)$0, and $x^{k+1}=x^k+f(x^k;\theta^k)$1, respectively; with Mapillary pre-training, SegNode + OCR reached $x^{k+1}=x^k+f(x^k;\theta^k)$2 mIoU on Cityscapes test (Khoshsirat et al., 2022).

In reduced-order modeling for environmental hydrodynamics, Neural ODEs were used to propagate latent dynamics in POD-based models for incompressible flow around a cylinder and shallow-water systems. Relative to POD + RBF and DMD non-intrusive reduced-order models, the Neural ODE formulation was reported as stable, accurate, and promising for extrapolatory prediction, though with a large training-time burden: cylinder runs required about $x^{k+1}=x^k+f(x^k;\theta^k)$3 h on a single V100 GPU, and more generally training ranged from $x^{k+1}=x^k+f(x^k;\theta^k)$4 k to $x^{k+1}=x^k+f(x^k;\theta^k)$5 k epochs (Dutta et al., 2021).

For irregular and noisy observations, an alternative formulation couples two networks: an implicit time-to-state representation $x^{k+1}=x^k+f(x^k;\theta^k)$6 and a vector-field network $x^{k+1}=x^k+f(x^k;\theta^k)$7. Training minimizes a weighted combination of data fit, an integral-form ODE consistency loss, and a gradient consistency loss. In a cubic oscillator experiment, the method substantially outperformed a standard Neural ODE baseline under increasing noise; for example, at $x^{k+1}=x^k+f(x^k;\theta^k)$8 noise the reported mean vector-field recovery error was $x^{k+1}=x^k+f(x^k;\theta^k)$9 for the implicit Neural ODE formulation versus $f(\cdot,t;\theta)$0 for the standard baseline (Goyal et al., 2022).

Hybrid continuous–discrete systems motivate event-aware generalizations. Neural Event ODEs replace a fixed terminal time by event times defined implicitly by a neural event function $f(\cdot,t;\theta)$1, interrupt the continuous flow at each event, apply a discrete update $f(\cdot,t;\theta)$2, and resume integration. The method was evaluated on switching dynamical systems, multi-body collision, temporal point processes, and discrete control; on a switching linear system, the reported MSE values were $f(\cdot,t;\theta)$3 for an LSTM, $f(\cdot,t;\theta)$4 for a Neural ODE, and $f(\cdot,t;\theta)$5 for a Neural Event ODE (Chen et al., 2020).

4. Expressivity, robustness, generalization, and geometric limits

Because Neural ODEs define flow maps of ODEs, they inherit structural constraints from dynamical systems. The invertibility of the ODE flow is central in continuous normalizing flows and in many theoretical analyses, but it also limits representational capacity in low-dimensional settings: the survey literature notes that one may need to augment the state, for example by padding with zeros, in order to represent non-invertible mappings (Ruthotto, 2024).

Several works sharpen this point into exact embeddability and non-embeddability results. For the basic Neural ODE architecture $f(\cdot,t;\theta)$6 with output $f(\cdot,t;\theta)$7, the induced map is a homeomorphism under mild continuity and uniqueness assumptions, implying in one dimension that any realizable map must be strictly increasing; consequently, $f(\cdot,t;\theta)$8 is not embeddable by a basic Neural ODE. Adding a final linear layer increases capability but does not eliminate all obstructions. Augmenting the phase space enlarges the set of realizable transformations, yet maps with scalar components that are topological Morse functions with critical points remain non-embeddable in the purely augmented architecture. By contrast, combining augmentation with a final linear layer yields a universal embedding result for any Lebesgue-integrable $f(\cdot,t;\theta)$9, using $1$0 dimensions, and a suspension-flow construction shows that any smooth diffeomorphism can be exactly realized in $1$1 dimensions with two deterministic layers around the ODE block (Kuehn et al., 2023).

Time dependence is one route to additional expressivity. A non-autonomous Neural ODE takes the form

$1$2

with explicitly time-varying weights. The time-dependence study states that Neural ODEs are universal approximators only when they are non-autonomous, and proposes non-parametric weight trajectories based on bucketed, polynomial, trigonometric, random-feature, and hypernetwork parameterizations. In practice, trigonometric polynomial parameterizations were emphasized because they provide a clear degree-$1$3 control over expressiveness and smoothness, and reported CIFAR results include $1$4 for T-NANODE-10 on CIFAR-10 and $1$5 on CIFAR-100, compared with $1$6 and $1$7 for the autonomous NODE baseline at comparable activation memory of about $1$8 GB (Davis et al., 2020).

Robustness is another property linked to flow geometry. A theoretical analysis of ODENets attributes improved resistance to perturbations to the non-intersection of integral curves: if $1$9 is continuous in time and Lipschitz in the state, trajectories initialized at different points do not intersect. Empirically, ODENets outperformed parameter-matched CNNs on Gaussian and adversarial perturbations. For example, when trained on clean images only, ODENet accuracy on MNIST under Gaussian noise with $\phi_\theta:x\mapsto h(1)$0 and FGSM-$\phi_\theta:x\mapsto h(1)$1 was $\phi_\theta:x\mapsto h(1)$2 and $\phi_\theta:x\mapsto h(1)$3, versus $\phi_\theta:x\mapsto h(1)$4 and $\phi_\theta:x\mapsto h(1)$5 for the CNN. The same work introduced the time-invariant steady Neural ODE (TisODE), which adds a steady-state regularizer and removes explicit time dependence; on MNIST with clean + Gaussian-noised training data, TisODE improved ODENet performance from $\phi_\theta:x\mapsto h(1)$6 to $\phi_\theta:x\mapsto h(1)$7 under PGD-$\phi_\theta:x\mapsto h(1)$8 while maintaining comparable clean accuracy (Yan et al., 2019).

Statistical learning theory for Neural ODEs reflects the same structural emphasis on parameter size and temporal smoothness. For a family of continuous-time parameterized ODEs with time-varying coefficients bounded by $\phi_\theta:x\mapsto h(1)$9 and Lipschitz-in-time constant $\phi_\theta$0, a Lipschitz-covering-number argument yields a generalization bound whose leading terms scale as $\phi_\theta$1 plus an $\phi_\theta$2 contribution arising from the infinite-dimensional Lipschitz parameter class. The corresponding ResNet specialization produces a depth-independent bound involving the magnitude of the difference between successive weight matrices, $\phi_\theta$3, and numerical experiments on depth-$\phi_\theta$4, width-$\phi_\theta$5 MNIST ResNets showed a clear positive correlation between the measured inter-layer difference and the generalization gap (Marion, 2023).

5. Structured priors, uncertainty quantification, and stiff dynamics

A major development after the original Neural ODE paper has been the explicit incorporation of mechanistic structure. One line of work derives Lie symmetries and conservation laws not only for the forward hidden-state dynamics but also for the adjoint equations. In a simple Neural ODE with $\phi_\theta$6, symmetry analysis yields one forward conservation law $\phi_\theta$7 and three backward conservation laws $\phi_\theta$8, $\phi_\theta$9, and $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$00. These are inserted into the objective as squared-residual penalties:

$$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$01

In the charged-particle sinusoidal-field example, the symmetry-regularized model recovered $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$02 more accurately than an unregularized Neural ODE, exhibited smaller test error over held-out initial conditions, and displayed improved step-to-step stability in the ODE solver (Hao, 2023).

A second line addresses uncertainty quantification. One study applies the Laplace approximation to Neural ODEs, placing a Gaussian prior on the weights and approximating the posterior by a Gaussian centered at the MAP estimate with covariance given by an inverse Hessian or generalized Gauss–Newton approximation. The work considers vanilla MLP vector fields, symplectic Hamiltonian Neural ODEs, and physics-augmented models of the form $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$03. Reported experiments show that uncertainty bands expand outside the training domain, that mechanistic structure can collapse uncertainty along conserved directions, and that uncertainty reveals the effect of mechanistic priors more directly than point predictions alone (Ott et al., 2023).

A broader Bayesian program treats the Neural ODE weights as full posterior random variables rather than as Gaussian approximations around a point estimate. The Bayesian Neural ODE study evaluates NUTS, SGHMC, SGLD, and variational inference with normalizing flows, and reports successful integration with classical physical systems, MNIST, and universal ODE system identification. On MNIST, the abstract reports a posterior sample accuracy of $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$04 on the test ensemble of $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$05 images (Dandekar et al., 2020).

Stiffness poses a separate but central challenge. For Robertson’s problem with time scales from $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$06 to $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$07 and for the $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$08-species POLLU mechanism, naïve use of a stiff solver and the continuous adjoint was found to fail unless three ingredients were combined: a deeper network with rectified activations such as GELU, scaling of both the ODE and the loss so that variables of very different magnitudes contribute comparably, and stabilized gradient computation through either a discrete adjoint or a quadrature-adjoint method that avoids integrating the stiff problem backward in a single shot. Under these modifications, a scaled $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$09-GELU Neural ODE on Robertson drove MAE below $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$10 within $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$11 epochs and achieved relative errors below $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$12 for all three species over five decades of time, while the POLLU model achieved MAE below $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$13 (Kim et al., 2021).

A complementary approach to stiffness is data-driven time reparametrization for model order reduction. Using adaptive time-stepping from an implicit Radau IIA reference solver, one defines a pseudo-time $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$14 that spaces trajectory samples uniformly over solver steps, learns a nonstiff ODE in $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$15, and separately learns the mapping between state and physical time. Across Van der Pol, OREGO, ROBER, E5, and POLLU benchmarks, the learned reparametrization reduced the cost of online inference, with the conclusion reporting $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$16–$$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$17 fewer right-hand-side calls and $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$18–$$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$19 faster wall-times while maintaining $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$20–$$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$21 accuracy and generalizing beyond the training interval (Caldana et al., 2024).

Structure can also be imposed directly for stability guarantees in discrete-time neural ODE analogues. A constrained formulation for discrete ordinary differential equations with algebraic nonlinearities enforces spectral-radius control of the learned transition matrix through row-wise softmax normalization and a damping mask, guaranteeing $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$22 by Perron–Frobenius arguments. Additional state and input inequality constraints are handled by ReLU-based slack penalties. On a bilinear thermal building benchmark, constrained models consistently outperformed unconstrained ones in long-horizon open-loop prediction; at horizon $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$23, the reported open-loop MSE was about $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$24 for the unconstrained black-box model and about $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$25 for the constrained black-box model (Tuor et al., 2020).

6. Solver-aware architectures, interpretability, and open limitations

Neural ODE performance depends not only on the learned vector field but also on the numerical method used to evaluate it. A meta-solver approach treats the solver itself as a trainable source of variability by sampling from parameterized Runge–Kutta families during training. For the two-stage second-order family, the Butcher tableau is parameterized by a scalar $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$26, interpolating between midpoint and Heun’s method. The study argues that training with a distribution over solvers reduces overfitting to a single discretization and can improve adversarial robustness. On CIFAR-10 at $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$27, standard training with RK2($$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$28) yielded FGSM and PGD robust accuracies of $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$29 and $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$30, while solver smoothing with Gaussian perturbations around $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$31 yielded $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$32 and $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$33 (Gusak et al., 2021).

Interpretability has recently been pushed beyond first-order sensitivity analysis. A 2025 framework introduces Event Transition Tensors (ETTs), high-order Taylor coefficients for event maps associated with Neural ODE trajectories reaching event manifolds. The method expands the event map

$$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$34

to order $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$35 around a nominal state-parameter pair and uses multivariate Faà di Bruno expansions to propagate uncertainty analytically rather than by Monte Carlo. Applications include a Neural ODE prey–predator control model, neural optimal feedback for comet landing, and a three-body neural Hamiltonian system, with the stated goal of expressing learned dynamics through explicit mathematical structures rather than opaque black-box behavior (Izzo et al., 2 Apr 2025).

At the same time, several limitations remain stable across the literature. Solver overhead and function-evaluation counts can exceed those of fixed-depth ResNets, particularly on stiff dynamics; the continuous formulation may under-parametrize mappings that are easier in the discrete case; adjoint-based training can suffer from numerical instability; and activation recomputation, state augmentation, and asynchronous events complicate batching and large-scale deployment (Ruthotto, 2024, Chen et al., 2020). Domain-specific studies make the same point in more concrete terms: semantic segmentation shows large memory savings but still requires tolerance tuning, environmental reduced-order modeling reports training times of $$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$36–$$\frac{d\mathbf{h}(t)}{dt}=f\bigl(\mathbf{h}(t),t;\theta\bigr),$$37 h/GPU per architecture, and stiff-system studies require specialized scaling and adjoint stabilization before training becomes reliable (Khoshsirat et al., 2022, Dutta et al., 2021, Kim et al., 2021).

The cumulative picture is therefore twofold. Neural ODEs furnish a mathematically coherent continuous-time generalization of residual networks, with adaptive computation, exact flow-based density modeling, and strong connections to dynamical systems, PDEs, and optimal control (Chen et al., 2018, Ruthotto, 2024). At the same time, their practical success depends on solver–model co-design, control of stiffness and reversibility, and architectural choices that reconcile flow geometry with target expressivity. This suggests that the enduring significance of Neural ODEs lies not in a single architecture, but in a family of continuous-depth methods whose properties are jointly determined by vector-field parameterization, numerical analysis, and the structure of the underlying scientific or statistical task.