Neural ODEs: Learning Continuous-Time Dynamics with Neural Networks

A Cup of Chai and a Question About Layers

Set a hot cup of chai on a table in an air-conditioned room and watch a thermometer in it. The temperature does not fall in jumps — it does not sit at 90°C for ten minutes and then instantly drop to 80°C. It falls continuously, a smooth curve bending toward room temperature. Newton's law of cooling describes this with a rate equation: the speed at which temperature drops is proportional to how much hotter the chai still is than the room. If h(t) is the temperature at time t, this law is written

dh/dt = -k * (h - h_room)

This is a differential equation: instead of telling you h(t) directly, it tells you the rate of change of h at every instant, and you must recover the full curve from that rate. Every simulation you have ever seen — a ball falling, a satellite orbiting, an epidemic spreading, a voltage decaying in a capacitor — is a version of this same idea: a rule for the instantaneous rate of change, from which the whole trajectory must be reconstructed.

Now compare that to a deep neural network. A network with 50 layers takes an input h_0 and transforms it 50 times, output of one layer feeding the next: h_1, h_2, h_3, …, h_50. That is not a continuous curve — it is a discrete staircase, 50 abrupt jumps. This chapter is about a genuinely elegant idea, published by Chen, Rubanova, Bettencourt, and Duvenaud in 2018: if you let the number of layers go to infinity and the jump at each layer shrink to zero, a deep network's staircase turns into a smooth curve governed by a differential equation — except now a neural network supplies the rate-of-change rule instead of a physicist. This is a Neural ODE: a model that does not learn a fixed number of transformations, but learns the continuous law of motion of a hidden state, and lets a numerical solver figure out how far to run it.

To understand this properly you need three pieces, in order: how to numerically solve a differential equation when you know the rate-of-change rule (Euler's method); why a residual network is secretly already doing exactly that; and how you can train a network whose "layers" are really steps of a solver, without storing every single step in memory (the adjoint sensitivity method). We will derive all three, not just describe them.

Solving an ODE One Small Step at a Time

Take the simplest possible rate law: dh/dt = -h, with h(0) = 1. This says "the state shrinks at a rate equal to its current size" — pure exponential decay, the same shape as radioactive decay or the chai cooling example with k=1 and room temperature 0. This particular equation is solvable exactly by separating variables, a technique from the CBSE Class 12 differential equations chapter:

dh/h = -dt
∫ dh/h = -∫ dt
ln(h) = -t + C
h(t) = h(0) * e^(-t) = e^(-t)

So the exact answer at t=1 is h(1) = e^(-1) = 0.3678794…. Knowing the exact answer is what lets us grade an approximate method honestly, which is exactly what we do next.

Most differential equations that show up in real machine learning problems are not solvable by a clean formula like this — especially once the rate-of-change rule is a neural network with thousands of tangled parameters instead of a simple -h. So instead of solving exactly, we approximate the curve by taking small, straight-line steps, using the one piece of information the differential equation actually gives us: the slope at the current point.

This is Euler's method. If you know h at time t, and you know the rate dh/dt = f(h,t) at that instant, then over a very short interval dt the curve barely bends, so you can walk forward along the tangent line:

h(t + dt) ≈ h(t) + dt * f(h(t), t)

This comes directly from the definition of a derivative: f(h,t) = dh/dt ≈ [h(t+dt) - h(t)] / dt, and multiplying both sides by dt and rearranging gives the update rule above. Repeat this update in a loop, re-evaluating the slope at each new point, and you trace out an approximate curve made of tiny straight segments — a staircase that hugs the true curve more and more tightly as the steps get smaller.

For our decay equation, f(h,t) = -h, so one Euler step is simply h ← h - dt*h = h*(1-dt). Starting from h=1 and taking n steps of size dt = 1/n to reach t=1 gives the closed form h(1) ≈ (1-1/n)^n. Here is that computation, traced exactly, for increasing numbers of steps:

def euler_decay(n_steps):
    h = 1.0
    dt = 1.0 / n_steps
    for _ in range(n_steps):
        h = h + dt * (-h)   # f(h, t) = -h
    return h

for n in [10, 20, 40, 100]:
    print(n, euler_decay(n))
# 10  0.3486784401000001
# 20  0.3584859224085419
# 40  0.36323243988788034
# 100 0.3660323412732292

Compare each result against the true value e^(-1) = 0.3678794…:

Two things are worth deriving rather than just observing. First, why is there error at all? Expand the true solution with a Taylor series around t:

h(t+dt) = h(t) + dt*h'(t) + (dt²/2)*h''(t) + …

Euler's method keeps only the first two terms and throws away everything from (dt²/2)*h''(t) onward. So the error introduced in a single step is proportional to dt². But reaching t=1 takes 1/dt steps, and these per-step errors accumulate, so the total error at the end scales like (1/dt) × dt² = dt — proportional to dt itself, not dt². This is why Euler's method is called a first-order method: global error shrinks in direct proportion to the step size.

Second, does the table actually show that? Look at the first three rows, where the step count doubles cleanly each time (10 → 20 → 40): the error goes 0.019201 → 0.009394 → 0.004647, each roughly half the one before — exactly what "error proportional to dt" predicts, since dt itself is halving. The 100-step row is not a clean doubling of the 40-step row (it is a further, smaller refinement), but it continues the same downward trend, landing closest to the true answer at a relative error of about 0.50%. This is the sense in which smaller steps make an ODE solver more faithful to the real continuous curve — a fact we will use constantly for the rest of the chapter, since "how many solver steps" is one of the actual design choices in a Neural ODE.

Where the Neural Network Enters

Euler's method needed one thing we simply assumed: the rate function f(h,t). For cooling chai, Newton gave us f = -k(h - h_room) from physics. But consider a harder problem: you have irregular blood-glucose readings from a patient taken every few hours, or you have a population's disease-spread numbers logged on inconsistent dates, or you have a satellite's tracked position at unevenly spaced ground-station contacts. You suspect the underlying process obeys some continuous rate law, but you have no physics textbook that hands you the formula.

A Neural ODE's core move is to stop assuming you know f, and instead let a neural network learn it: replace the rate function with f(h,t,θ), a neural network with trainable weights θ, so that

dh/dt = f(h, t, θ)

To predict, you pick a starting hidden state h(t₀) (say, an encoding of the patient's first reading), then solve this ODE forward in time using any numerical solver — Euler's method, or something more accurate — to get h(t₁), the state at whatever later time you care about. To train, you compare that predicted state to the real observation and adjust θ so the network's learned dynamics better match reality. The solver can be asked for the state at any time, not just fixed integer steps — which is precisely why Neural ODEs handle irregularly-timed data naturally, unlike a standard recurrent network, which is built around one fixed step between readings and struggles when the gaps between real observations vary.

The Reveal: A ResNet Layer Is One Euler Step

Here is the connection that makes Neural ODEs feel less like a new invention and more like an inevitable discovery. A residual network (ResNet), one of the most successful deep architectures, stacks layers with a "skip connection": instead of the next layer fully replacing the hidden state, it only adds a correction to it.

h_(t+1) = h_t + F(h_t, θ_t)

Now look again at the Euler update we derived two sections ago:

h(t + dt) = h(t) + dt * f(h(t), t, θ)

Set dt = 1 and rename F = f. The two equations are identical. A ResNet layer is one Euler step of an ODE solver, with step size fixed at exactly 1, applied to a rate function that happens to be a different neural network at every layer. Stacking 50 ResNet layers is running Euler's method for 50 steps of size 1 through 50 different learned functions f_1, f_2, …, f_50.

Chen et al.'s question was: what if, instead of 50 different functions taking 50 steps of size 1, we used a single function f(h,t,θ) that also takes time t as an input, and let a solver take as many steps of whatever size it needs to integrate accurately from t=0 to t=1? You are no longer choosing a discrete "depth" (number of layers) at all — depth becomes a continuous integration variable, and how finely to divide it becomes the solver's job, not an architecture decision. This is the "continuous-time" and "continuous-depth" language used for Neural ODEs: time in the ODE literally plays the role that layer-index plays in a normal deep network.

The diagram below makes this concrete for our decay example. The smooth curve is the true continuous solution h(t) = e^(-t). The staircase is the 10-step Euler approximation from the table above — and each corner of that staircase is exactly what a "layer output" would be in a 10-layer residual network implementing this same dynamics.

As the number of layers grows and the step size shrinks, the red staircase hugs the blue curve ever more closely — which is exactly the convergence behaviour proven in the table. A Neural ODE is the limit of this process: infinitely many, infinitesimally small residual layers, replaced by one call to a numerical ODE solver.

Solving the Forward Pass: Which Solver?

Euler's method is the simplest possible solver and the easiest to derive by hand, which is why we used it above, but real Neural ODE implementations rarely use it in production because a first-order method (error proportional to dt) needs very small steps to be accurate, which means many function evaluations. Higher-order solvers such as the classical fourth-order Runge–Kutta method (RK4) use several slope evaluations per step, combined in a weighted average, to cancel out more terms of the Taylor expansion — RK4's global error is proportional to dt⁴, so doubling the step size only multiplies its error by roughly 16, not 2, letting it take far fewer, larger steps for the same accuracy. Better still, adaptive solvers (such as Dormand–Prince, the default in the original Neural ODE paper's implementation) automatically shrink dt in regions where f changes rapidly and grow it where the trajectory is nearly straight, so the "number of layers" a Neural ODE effectively uses is not fixed in advance at all — it is decided at runtime by how complicated the learned dynamics turn out to be. This is a genuine qualitative difference from an ordinary neural network, whose depth is a hard architectural choice made before training even starts.

Training Backward Through a Continuous Solver: The Adjoint Method

Training any neural network requires computing how the loss L depends on the parameters θ, via backpropagation. For an ordinary 50-layer ResNet, this means storing every one of the 50 intermediate hidden states during the forward pass, then walking backward layer by layer applying the chain rule. If a Neural ODE's solver takes, say, 500 adaptive steps to integrate accurately, naively backpropagating "through the solver" the same way would mean storing all 500 intermediate states — memory cost growing directly with solver accuracy, which is a bad trade since we specifically wanted the freedom to take as many steps as needed.

The 2018 paper's second key contribution is a way to avoid this entirely, called the adjoint sensitivity method. It is worth deriving properly, because the derivation also fixes an important shift we need to flag explicitly: up to now, the decay example treated h as a single number for clarity. In a real network, h is a hidden-state vector — a list of many numbers (say 64 or 256 hidden units) — and f(h,t,θ) is a vector-valued function that outputs one rate-of-change number per hidden unit. Because of this, the derivative of f with respect to h is no longer a single slope; it is a matrix of partial derivatives called the Jacobian, written ∂f/∂h, whose entry in row i, column j tells you how the rate of change of hidden unit i responds to a small change in hidden unit j. Everywhere below, multiplying a vector by ∂f/∂h means an ordinary matrix-vector product, and I denotes the identity matrix (the matrix version of "multiply by 1", with 1s down the diagonal and 0s elsewhere) — it appears because "no change" in the residual/Euler update below must be represented the same way for a whole vector as the number 1 represents it for a scalar.

Define the adjoint a(t) = ∂L/∂h(t): a vector, the same length as h(t), recording how sensitive the final loss is to the hidden state at time t. We want to find how a(t) evolves. Start from the discrete Euler picture, one step at a time, exactly as in a ResNet's backpropagation:

h_(n+1) = h_n + dt * f(h_n, t_n, θ)

Differentiating this update with respect to h_n gives its Jacobian:

∂h_(n+1)/∂h_n = I + dt * (∂f/∂h_n)

By the chain rule, the adjoint at step n is obtained from the adjoint at step n+1 by multiplying through this Jacobian — precisely how gradients flow backward through any layer:

a_n = a_(n+1) * (I + dt * ∂f/∂h_n)
    = a_(n+1) + dt * a_(n+1) * (∂f/∂h_n)

Rearrange to isolate the change in the adjoint across one step:

(a_(n+1) - a_n) / dt = -a_(n+1) * (∂f/∂h_n)

The left-hand side is precisely the difference quotient that defines a derivative. As dt → 0 (infinitely many, infinitesimally small layers — the Neural ODE limit), it becomes da/dt, giving the adjoint ODE:

da/dt = -a(t) * (∂f/∂h(t))

This says: the sensitivity of the loss to the hidden state obeys its own differential equation, running backward in time from t₁ (where a(t₁) = ∂L/∂h(t₁) is just the ordinary loss gradient at the output) down to t₀. To train the network, you solve this adjoint ODE backward with the same kind of numerical solver used for the forward pass, and along the way you accumulate the gradient with respect to the parameters θ using a closely related formula, dL/dθ = -∫ a(t)·(∂f/∂θ) dt, integrated over the same backward sweep.

The crucial practical payoff: to run this backward ODE, the solver needs to know h(t) at each point it visits — but it does not need those values pre-stored, because h(t) obeys its own ODE too, and can simply be re-solved backward in time alongside the adjoint, starting from the already-known h(t₁). So training needs only O(1) memory relative to the number of solver steps, regardless of how many steps the forward solve took — a sharp contrast with a standard deep network, whose training memory grows linearly with depth because every layer's activations must be kept around for the backward pass.

Where This Genuinely Gets Used

Two properties make Neural ODEs a real, distinct tool rather than a mathematical curiosity. First, irregular timing: because the solver can be queried for h at any real-valued time, Neural ODEs are a natural fit for medical time series where lab readings arrive at uneven intervals, or for satellite telemetry where ground-station contact windows are irregular and ISRO-style orbit-propagation models must interpolate a continuous trajectory between sparse tracking fixes. A standard RNN, built around one fixed step per token, has no clean way to represent "the next reading is 11 hours away" versus "3 hours away" without arbitrary padding or resampling; a Neural ODE just integrates for the right amount of time. Second, memory-efficient depth: because training cost does not scale with the number of solver steps, Neural ODEs let you trade compute for accuracy at inference time (take more steps, get a more faithful trajectory) without redesigning or retraining the network — something impossible for a network whose "depth" is a fixed count of learned layers baked in during training.

A Misconception Worth Correcting

A common misreading of "Neural ODE" is to assume it means the differential equation itself gets solved analytically, the way you solved dh/dt = -h by separating variables earlier in this chapter, and that the neural network is somehow a closed-form formula you could write out and integrate by hand. This is false, and it is precisely why the machinery of this chapter (Euler's method, adaptive solvers, the adjoint method) exists at all: once f(h,t,θ) is a neural network with nonlinear activation functions and thousands of parameters, the ODE almost never has a clean analytic solution. It is solved the same way you would numerically integrate any equation you cannot solve by hand — by stepping forward with a solver, exactly as in the table above, just with a far more complicated f. The "neural" in Neural ODE describes what supplies the rate-of-change rule, not a shortcut around numerical integration.

CBSE and Competitive-Exam Connections

The separation-of-variables solution of dh/dt = -h is a standard technique from the Class 12 differential equations unit, and identifying order (first) and degree (one, since the highest derivative appears to the first power) of the governing equation is directly examinable. The Taylor-series argument used to derive Euler's O(dt) global error is the same expansion technique tested in JEE/BITSAT calculus questions on approximation and error bounds. Numerical solving of differential equations by small steps (Euler's method and its refinements) appears explicitly in GATE-level numerical methods and is a natural KVPY/Olympiad-style question: "if a first-order method has global error proportional to dt, how must you scale the step count to halve your error?" — a question this chapter's convergence table answers directly (halve dt, i.e. double the steps). Finally, implementing the exact euler_decay function above is a legitimate CBSE Computer Science (Python) exercise: it only uses a loop, a running variable, and arithmetic, yet it is genuinely simulating a continuous physical process.

Active Recall

1. Solve dh/dt = -h by separation of variables yourself, starting from ∫dh/h = -∫dt, and confirm you reach h(t) = h(0)e^(-t).
2. Take dh/dt = 2h, h(0) = 1. Using Euler's method with dt = 0.25, compute h after 2 steps (i.e. at t = 0.5) by hand. Then compare against the exact solution h(t) = e^(2t) at t=0.5, and state the relative error as a percentage.
3. Explain, in one sentence each: (a) why a 50-layer ResNet is mathematically one specific choice of ODE solver and step size; (b) what changes if you instead solve the same underlying ODE with 500 adaptive steps instead of 50 fixed steps.
4. In the adjoint ODE derivation, why does the identity matrix I appear in ∂h_(n+1)/∂h_n = I + dt·∂f/∂h_n rather than just the number 1? What would this equation look like if h were a single scalar instead of a vector?
5. A classmate says, "Neural ODEs need more memory to train than a ResNet of equivalent depth, because they take hundreds of solver steps." Identify exactly what is wrong with this claim and name the technique that makes it false.

Summary

A differential equation specifies a rate of change, dh/dt = f(h,t), and Euler's method reconstructs the underlying curve by repeatedly stepping along the current tangent line, an approximation whose global error is provably proportional to the step size dt — verified directly in this chapter's convergence table, where halving the step count from 10 to 20 to 40 roughly halved the error each time, and taking it further to 100 steps brought the relative error down to about 0.50%. A residual network's layer update, h_(t+1) = h_t + F(h_t,θ_t), is exactly one Euler step with dt=1; letting the number of such steps go to infinity while their size shrinks to zero turns a stack of discrete layers into a continuous-time trajectory governed by dh/dt = f(h,t,θ), where f is now a neural network — this is a Neural ODE. Its forward pass is computed by any ODE solver, from simple Euler to adaptive Runge–Kutta variants; its backward pass, rather than storing every solver step, runs a second differential equation — the adjoint ODE, da/dt = -a(t)·(∂f/∂h), where a and h are vectors and ∂f/∂h is their Jacobian — backward in time, recomputing h(t) along the way instead of caching it, which is what gives Neural ODEs their signature O(1) training memory and their natural ability to handle irregularly timed data, from medical readings to satellite tracking fixes.