Interactive Differential Equations

Definition. An ordinary differential equation (ODE) is a relation between an unknown function \(x(t)\) and its derivatives. The most common first-order form is \(\frac{dx}{dt}=f(t,x)\), and the degree of the highest derivative is called the order of the equation. "Ordinary" means there is only one independent variable; with several it becomes a partial differential equation.

An equation, but the answer is a function. The answer to \(x^2-4=0\) is the number \(\pm 2\), but the answer to \(\dot{x}=2x\) is the function \(x(t)=Ce^{2t}\). Moreover there is one for every constant \(C\) — not a single solution but a family of functions. That is why several curves are drawn at once in the demo, and picking one out requires the initial condition of the next section.

A local rule builds a global trajectory. \(f(t,x)\) tells us only the slope at each point. Drawing those slopes densely as short segments on the plane gives the direction field, and a solution curve is a curve tangent to those segments at every instant. The important point is that without knowing how to solve the equation, the shape of its solutions can be read off this picture alone.

Where do they come from? Newton's second law \(m\ddot{x}=F\) is a second-order ODE, and the population model \(\dot{P}=kP\) is a direct transcription of "the rate of increase is proportional to the current population." The moment the phrase "rate of change" appears in physics, chemistry, biology or economics, a differential equation is being written.

In machine learning. In the limit of the learning rate going to 0, gradient descent becomes the ODE \(\dot\theta=-\nabla L(\theta)\) (gradient flow). A Neural ODE takes the network itself as \(f\), solves \(\dot{h}=f_\theta(h,t)\), and uses the solution in place of layers. And the generation process of the diffusion models this article is heading toward is, in the end, solving a single learned ODE.

\[\begin{aligned} \frac{dx}{dt} &= f(t, x) \\ \dot{x} = 2x &\;\Rightarrow\; x(t) = C e^{2t} \\ m\ddot{x} &= F, \qquad \dot{P} = kP \end{aligned}\]

Change the equation with the presets. The short segments in the background are the direction field, and the coloured curves are solutions corresponding to different constants C. Check that every solution curve is exactly tangent to the segment at its position — that is all the equation says.

Definition. An initial value problem (IVP) gives a differential equation together with an initial condition: \(\dot{x}=f(t,x),\; x(t_0)=x_0\). For an \(n\)-th order equation, \(n\) conditions are needed, from \(x(t_0)\) through \(x^{(n-1)}(t_0)\).

Geometrically: picking one curve. The direction field is densely paved with non-intersecting solution curves. An initial condition \((t_0,x_0)\) marks a single point on the plane, and the curve through that point is the unique solution. Dragging the point in the demo makes the solution curve slide along with it.

Existence and uniqueness. It does not work for any \(f\) whatsoever. The Picard–Lindelöf theorem says that if \(f\) is continuous and satisfies a Lipschitz condition in \(x\), then locally a solution exists and is unique. If the Lipschitz condition breaks, uniqueness can fail — \(\dot{x}=\sqrt{|x|},\,x(0)=0\) has both \(x\equiv 0\) and \(x=t^2/4\) as solutions. The presets let you see this case directly.

Finite-time blow-up. Even when a solution exists, it may not exist forever. The solution of \(\dot{x}=x^2,\,x(0)=1\) is \(x=1/(1-t)\), diverging to infinity as \(t\to 1\). The equation is perfectly well behaved while the solution vanishes in finite time, and this is one of the reasons values become NaN in numerical simulation.

In machine learning. In a Neural ODE the initial value is the input and the value at the final time is the output. Because uniqueness of the solution is guaranteed, this correspondence is well defined as a function, and thanks to the property that trajectories cannot cross, the map is automatically invertible — which is exactly the grounds on which a continuous-time normalizing flow holds.

\[\begin{aligned} \dot{x} &= f(t,x), \qquad x(t_0) = x_0 \\ \text{Lipschitz: } &|f(t,x)-f(t,y)| \le L|x-y| \\ &\Rightarrow \text{the solution is locally unique} \end{aligned}\]

Drag the initial point and one solution curve through it is selected. On the "uniqueness failure" preset, two different solutions extend from the same initial point, and on "finite-time blow-up" the solution shoots up vertically at a particular time.

Definition. If \(t\) does not appear explicitly on the right-hand side, as in \(\dot{x}=f(x)\), the system is autonomous; if it does, as in \(\dot{x}=f(t,x)\), it is non-autonomous. An autonomous system is the situation "the laws of physics are the same yesterday and today," while a non-autonomous one has a time-varying force coming in from outside.

An autonomous system is invariant under time translation. If \(x(t)\) is a solution then so is \(x(t+c)\). So "when did we set out" does not matter, only "where did we set out from." Thanks to this we may erase the time axis and draw only the state space with no loss of information — and that picture is the phase portrait of the next section.

Trajectories do not cross. In an autonomous system there is only one arrow at each point, so if two trajectories meet at a point they must coincide from then on. Hence distinct trajectories never cross. In a non-autonomous system, by contrast, the arrow at the same location differs by time, so drawing them in \(x\)-space alone rather than the \(x\)-\(t\) plane makes the trajectories appear to cross.

Turning a non-autonomous system autonomous. Promote time to part of the state. Adding a variable \(s=t\) and writing \(\dot{s}=1,\;\dot{x}=f(s,x)\) removes \(t\) from the right-hand side and gives an autonomous system. It is a trade: one extra dimension in exchange for a simpler theory.

In machine learning: why time is fed as an input. The noise-prediction network of a diffusion model always takes the time \(t\) as an input, as in \(\epsilon_\theta(x, t)\). That is because the forward noise schedule is a non-autonomous system varying with time, and even for the same \(x\) the job to be done differs entirely between an early and a late stage. Whether or not to feed \(t\) into \(f_\theta(h,t)\) in a Neural ODE bears directly on expressiveness for the same reason.

\[\begin{aligned} \dot{x} &= f(x) &&\text{autonomous} \\ \dot{x} &= f(t,x) &&\text{non-autonomous} \\ \dot{s} = 1,\; \dot{x} &= f(s,x) &&\text{made autonomous} \end{aligned}\]

Alternate between the two systems. In the autonomous system the slope at a given height is the same whatever the time, so each horizontal row of the direction field is identical. In the non-autonomous system the slope at the same height differs by time, so the shape of a solution curve depends on when it started.

Definition. For a two-dimensional autonomous system \(\dot{x}=f(x,y),\;\dot{y}=g(x,y)\), the picture drawing the vector field \((f,g)\) and representative trajectories on the state space \((x,y)\) is the phase portrait. The vector field of the previous article reappears here unchanged — except that its streamlines now carry the interpretation of "motion in time."

Fixed points. At points where \(f=g=0\) the arrows vanish. A state placed there never moves, so it is called a fixed point (an equilibrium). The long-run behaviour of a system is largely determined by "which fixed point does it get pulled toward."

Linearization gives the classification. Approximating the system near a fixed point by its Jacobian gives \(\dot{\delta}\approx J\delta\), and the eigenvalues of \(J\) determine the type. Two negative real eigenvalues give a stable node (drawn in), two positive ones an unstable node, and mixed signs a saddle (attracted along one direction and repelled along another). Complex eigenvalues mix in rotation to give a spiral, with the sign of the real part deciding whether it winds in or unwinds. A single eigenvalue from linear algebra reads off the fate of the dynamics.

Why the eigenvalues decide the fate. The solution of \(\dot\delta=J\delta\) is of the form \(\delta(t)=e^{\lambda t}v\) along an eigenvector direction. That is, along an eigenvector the motion reduces to a single exponential. So if \(\lambda<0\) that component dies exponentially (drawn in), and if \(\lambda>0\) it grows exponentially (pushed away). A general starting point is a sum of eigenvectors, so each direction shrinks or grows separately according to its own \(\lambda\), and ultimately the direction with the largest (least negative) \(\lambda\) dominates the picture in the long run. A trajectory near a saddle approaching along one direction and eventually shooting out along another is the result of this competition.

Why complex eigenvalues mean rotation. For \(\lambda=a\pm bi\) we have \(e^{\lambda t}=e^{at}(\cos bt \pm i\sin bt)\), so the real part \(a\) controls the magnitude while the imaginary part \(b\) rotates the angle through \(\cos\) and \(\sin\). The two roles multiply to give a spiral that "turns while growing or shrinking." The condition for rotation can be read geometrically too — having a real eigenvector means there exists a line whose direction the transformation preserves, whereas a pure rotation has no such line at all. So the case with no real eigenvector (complex eigenvalues) is exactly the case with rotation mixed in. When \(a=0\) the magnitude does not change and the orbit is a closed ellipse, which is precisely the situation in which GAN training fails to converge and oscillates instead.

Limit cycles. Fixed points are not the only destinations. Some systems converge to a closed orbit, called a limit cycle. Self-sustaining rhythms such as a heartbeat or predator–prey oscillation have this structure.

In machine learning. The fixed points of gradient flow \(\dot\theta=-\nabla L\) are the critical points, and its Jacobian is the Hessian. A positive definite Hessian gives a stable node (a local minimum) and mixed signs a saddle point — the critical-point test of the previous article rewritten in the language of dynamics. Meanwhile a system like a GAN, where two objectives compete, has a non-symmetric Jacobian with complex eigenvalues, and so shows spiral oscillation instead of convergence.

\[\begin{aligned} \dot{x} &= f(x,y), \quad \dot{y} = g(x,y) \\ \dot{\delta} &\approx J\delta \quad \text{(near a fixed point)} \\ \lambda_{1,2} &\;\Rightarrow\; \text{node · saddle · spiral} \end{aligned}\]

Click the figure to draw the trajectory starting from that point. The character of the fixed point differs by preset, and the readout shows the eigenvalues of the Jacobian and the type determined from them. On the predator–prey preset the trajectory draws a closed loop and circulates forever.

Definition. \(x_{n+1}=x_n+\Delta t\, f(t_n,x_n)\). Read the slope at the current position, travel straight in that direction for \(\Delta t\), and read the slope again where you land. Repeat.

It is a truncated first-order Taylor expansion. Keeping only up to the first-order term of \(x(t+\Delta t)=x(t)+\Delta t\,\dot{x}(t)+\frac{\Delta t^2}{2}\ddot{x}+\cdots\) gives the Euler method. The Taylor expansion of the previous article becomes an algorithm here — and what was discarded is the error.

Local and global error. The error of one step is \(O(\Delta t^2)\), but covering the interval \([0,T]\) takes \(T/\Delta t\) steps, so the accumulated error is \(O(\Delta t)\). Hence the Euler method is called first-order accurate. Halving \(\Delta t\) roughly halves the error, and this slow convergence is what motivated the search for better methods.

Geometrically, why it drifts off. The Euler method uses only the slope at the start of each step. But the slope changes within the interval, so it is always pushed outward on the side toward which the curve bends. Increasing \(\Delta t\) in the demo shows the approximate trajectory departing systematically to the outside of the true solution — a bias, not a random error.

In machine learning. The deterministic samplers of diffusion models, DDIM among them, are essentially the Euler method (or a variant) applied to the probability flow ODE. Reducing the number of sampling steps means increasing \(\Delta t\), and this truncation error then shows up directly as degraded image quality. The reason higher-order samplers such as DPM-Solver appeared is in the next section.

\[\begin{aligned} x_{n+1} &= x_n + \Delta t\, f(t_n, x_n) \\ \text{local error} &= O(\Delta t^2) \\ \text{global error} &= O(\Delta t) \quad \text{(first-order accuracy)} \end{aligned}\]

Set Δt large and the stepwise approximation (orange) departs noticeably from the true solution (grey). The short tangent drawn at each step shows the rule "read the slope at this point and go straight" exactly. Check in the readout whether halving Δt roughly halves the error.

Definition (RK4). Each step measures the slope at four places and takes a weighted average: \(k_1\) at the start, \(k_2\) at the midpoint estimated by \(k_1\), \(k_3\) at the midpoint re-estimated by \(k_2\), and \(k_4\) at the end. The final update is \(x_{n+1}=x_n+\frac{\Delta t}{6}(k_1+2k_2+2k_3+k_4)\).

Why the midpoint twice? The flaw of the Euler method was representing a whole interval by the slope at its start. Using the slope in the middle of the interval gives a much better representative (the midpoint method, second order), and averaging several points with weights cancels still higher terms of the Taylor expansion. The weights \((1,2,2,1)/6\) are the same coefficients as Simpson's rule for integration — not a coincidence but the same principle.

The power of order. RK4 has global error \(O(\Delta t^4)\). Halving \(\Delta t\) divides the error by sixteen. Each step costs four function evaluations, so a step is four times as expensive, but the number of steps needed for the same accuracy is incomparably smaller, so in practice it is far cheaper. Comparing the three methods at the same \(\Delta t\) in the demo opens this difference by orders of magnitude.

It is not free. Higher-order methods pay off only when the solution is smooth. If \(f\) changes abruptly or is discontinuous, the assumptions behind the high-order terms break and the order collapses. Also, for the stiff problems of the next section, \(\Delta t\) cannot be increased however high the order — accuracy and stability are separate matters.

In machine learning. Neural ODEs internally use adaptive RK-family solvers (Dormand–Prince and others), so the "number of layers" is not fixed but determined by the difficulty of the solution. Fast samplers for diffusion models such as DPM-Solver and DEIS are also higher-order methods exploiting the special structure of the probability flow ODE (integrating the linear term exactly and approximating only the nonlinear part), which is how they cut 1000 steps down to 10–20.

\[\begin{aligned} k_1 &= f(t_n, x_n) \\ k_2 &= f(t_n + \tfrac{\Delta t}{2}, x_n + \tfrac{\Delta t}{2}k_1) \\ k_3 &= f(t_n + \tfrac{\Delta t}{2}, x_n + \tfrac{\Delta t}{2}k_2) \\ k_4 &= f(t_n + \Delta t, x_n + \Delta t\,k_3) \\ x_{n+1} &= x_n + \tfrac{\Delta t}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{aligned}\]

Compare the three methods at the same Δt. In the log-scale error graph below, the slope is the order of accuracy — Euler descends with slope 1, the midpoint method with 2, and RK4 with 4. Check that halving Δt divides the RK4 error by sixteen.

The test equation. The true solution of \(\dot{x}=\lambda x\;(\lambda<0)\) is \(x(t)=x_0e^{\lambda t}\), converging to 0. Applying the Euler method gives \(x_{n+1}=(1+\lambda\Delta t)x_n\), so \(x_n=(1+\lambda\Delta t)^n x_0\). For this to go to 0 we need \(|1+\lambda\Delta t|<1\), that is \(\Delta t < 2/|\lambda|\).

What happens past the boundary? The moment \(\Delta t\) exceeds \(2/|\lambda|\), the absolute value of the multiplier passes 1 and the sign flips every step while the magnitude grows — the true solution decreases monotonically while the numerical one bounces in alternating signs and explodes. Pushing \(\Delta t\) past the critical value in the demo produces this immediately. No amount of extra accuracy solves it.

Stiffness. Trouble arises when a system contains both a very fast and a very slow component (a large spread in eigenvalue magnitudes). What we care about is the slow behaviour, yet the step size is bound by the stability limit imposed by the fastest component. Such problems are called stiff, and explicit methods are effectively unable to solve them.

Implicit methods. Backward Euler, \(x_{n+1}=x_n+\Delta t\,f(t_{n+1},x_{n+1})\), uses the slope at the destination. For the test equation the multiplier becomes \(1/(1-\lambda\Delta t)\), whose absolute value is below 1 for \(\lambda<0\) however large \(\Delta t\) is (unconditionally stable). In exchange an equation must be solved at every step, which costs — stability bought with computation.

In machine learning. The upper bound on the learning rate in gradient descent comes from exactly the same calculation. For a quadratic loss, \(\theta_{n+1}=(I-\eta H)\theta_n\), so the stability condition is \(\eta<2/\lambda_{\max}(H)\). The rule of thumb "if the learning rate is too large the loss oscillates and diverges" is that inequality itself. And the condition number of the previous article reappears here — \(\lambda_{\max}\) binds the learning rate while \(\lambda_{\min}\) sets the speed of convergence, so a large condition number means a stiff problem.

\[\begin{aligned} \text{forward Euler: } &|1+\lambda\Delta t| < 1 \;\Rightarrow\; \Delta t < 2/|\lambda| \\ \text{backward Euler: } &\left|\tfrac{1}{1-\lambda\Delta t}\right| < 1 \;\; \forall \Delta t > 0 \\ \text{gradient descent: } &\eta < 2/\lambda_{\max}(H) \end{aligned}\]

Push Δt past the critical value (the red line). The true solution (grey) heads quietly to 0 while forward Euler (orange) bounces in alternating signs and explodes. Backward Euler (violet) is fine at the same Δt. Making λ more negative narrows the critical value — that is stiffness.

Definition. The flow of an autonomous system \(\dot{x}=f(x)\) is the map \(\Phi_t: x_0 \mapsto x(t)\) carrying an initial state to the state at time \(t\). The shift in perspective is regarding it as a transformation moving the whole space at once, rather than a single trajectory.

Group structure. A flow satisfies \(\Phi_0=\mathrm{id}\) and \(\Phi_s\circ\Phi_t=\Phi_{s+t}\) — flowing for \(t\) and then for a further \(s\) is the same as flowing for \(s+t\). Hence \(\Phi_t^{-1}=\Phi_{-t}\): every flow can be undone by running time backwards. This invertibility is used repeatedly throughout this article.

How does volume change? The volume of a region carried along by the flow is governed by the divergence: \(\frac{d}{dt}\log V=\nabla\cdot f\). A vector field with zero divergence (incompressible) preserves volume exactly, a positive one expands it, and a negative one contracts it. You can see this in the demo as a square is deformed along the flow.

In machine learning: Neural ODEs. Taking \(f\) to be a neural network makes \(\Phi_1\) behave like a single layer. The starting point was the observation that ResNet's \(h_{n+1}=h_n+f_\theta(h_n)\) is exactly one Euler step with \(\Delta t=1\) — a residual connection was a discretized flow all along. Instead of a layer count, the integration interval and the solver tolerance set the computation, and memory becomes independent of depth (the adjoint method).

Continuous normalizing flows. Because a flow is invertible, we can flow a random variable to transform its distribution. The log-density then follows \(\frac{d\log p}{dt}=-\nabla\cdot f\), and the decisive point is that only a trace need be computed rather than a Jacobian determinant. The Probability Flow ODE of the final section has exactly this structure.

\[\begin{aligned} \Phi_t &: x_0 \mapsto x(t) \\ \Phi_s\circ\Phi_t &= \Phi_{s+t}, \qquad \Phi_t^{-1} = \Phi_{-t} \\ \frac{d}{dt}\log V &= \nabla\cdot f, \qquad \frac{d\log p}{dt} = -\nabla\cdot f \end{aligned}\]

Move the time t and the whole grid, along with the square inside it, deforms along the flow. In a divergence-free field the square distorts while keeping its area, and in a contracting field the area shrinks. Turn t negative and it returns exactly to its original place — the invertibility of a flow.

Definition. A partial differential equation (PDE) is a relation between an unknown function of several variables and its partial derivatives. Typically one writes, for \(u(t,x)\), something like \(\frac{\partial u}{\partial t}=\alpha\frac{\partial^2 u}{\partial x^2}\) (the heat equation).

How does it differ from an ODE? The solution of an ODE is a single curve, while the solution of a PDE is an entire function deforming over time. If an ODE is "the motion of one point," a PDE is "the evolution of a whole shape." So the initial condition is not one value but an entire initial function \(u(0,x)\), and if space is finite, boundary conditions are needed too.

Three archetypes. Second-order linear PDEs fall broadly into three classes. Parabolic (the heat equation \(u_t=\alpha u_{xx}\)) smooths, loses information, and is hard to run backwards in time. Hyperbolic (the wave equation \(u_{tt}=c^2u_{xx}\), the transport equation) propagates at finite speed while preserving shape, and is reversible. Elliptic (Laplace's equation \(\Delta u=0\)) describes a timeless equilibrium. What we will meet in this article are the first two and their combination.

Mostly unsolvable by hand. Far more so than ODEs. So in practice we chop space into a grid, replace partial derivatives by differences, and integrate the resulting enormous ODE system in time (the method of lines). The demos are computed in exactly this way.

In machine learning. The protagonist of this article, the Fokker–Planck equation, is a PDE, and it describes how a probability density deforms over time. Understanding diffusion models is, in the end, understanding the relationship between this PDE and the SDE that produces it. Meanwhile the PINN and Neural Operator families, which solve physical equations themselves with neural networks, are an active area.

\[\begin{aligned} \frac{\partial u}{\partial t} &= \alpha\frac{\partial^2 u}{\partial x^2} &&\text{parabolic (diffusion)} \\ \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} &= 0 &&\text{hyperbolic (transport)} \\ \Delta u &= 0 &&\text{elliptic (equilibrium)} \end{aligned}\]

Apply the three equations to the same initial shape and compare. Transport carries the shape unchanged, diffusion blurs it, and the wave equation splits it left and right. Rewinding the time slider reveals that only diffusion has lost information and cannot restore the original shape.

Definition. \(\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=0\). For a constant speed \(c\) the solution is \(u(t,x)=u_0(x-ct)\), the initial shape pushed to the right by \(ct\). The shape does not change in the slightest.

Characteristics. Why this is so becomes clear by riding along with the flow. Moving along \(x(t)=x_0+ct\) and watching \(u\) gives \(\frac{d}{dt}u(t,x(t))=u_t+c\,u_x=0\), so to that observer \(u\) is constant. Such curves are called characteristics, and this structure — a PDE reducing to an ODE along characteristics — is the heart of hyperbolic equations.

When the speed varies with position. If \(c=c(x)\), the characteristics bend and the shape is stretched or compressed. When characteristics collide (a fast part behind catching a slow part in front) the solution begins to break (a shock), which is the difficulty of nonlinear conservation laws.

It loses no information. Transport is completely reversible. Flip the sign of \(c\), flow for the same time, and you return exactly to the original shape. This is the decisive property contrasting with the diffusion of the next section, and this difference returns in the final section as the question "is deterministic generation possible?"

In machine learning. The transport equation is a special case of the continuity equation of the next section (density conserved, velocity field constant). And what Flow Matching does is precisely a generalization of it — learning the velocity field carrying the data distribution to a target distribution. It is also where optimal transport theory meets generative models.

\[\begin{aligned} \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} &= 0 \\ u(t,x) &= u_0(x - ct) \\ \frac{dx}{dt} = c &\;\Rightarrow\; u \text{ is constant along that curve} \end{aligned}\]

Change the speed c and play. In the space–time picture below, the slanted stripes are the characteristics and their slope is the speed. Negative c flows leftward, and rewinding time restores the shape perfectly — transport loses no information.

Definition. For a density \(\rho(t,x)\) and a velocity field \(v(t,x)\), \(\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho v)=0\). Here \(j=\rho v\) is called the flux, and the equation reads "the time change of density + the divergence of the flux = 0."

Why divergence? Consider a small box. The only way the total inside can fall is by flowing out through the boundary. The net outflow through the boundary was exactly the divergence \(\nabla\cdot j\), so conservation becomes \(\partial_t\rho=-\nabla\cdot j\). The "sources and sinks" picture of the previous article acquires its physical meaning here.

The total is always conserved. Integrating over all space makes the divergence term 0 (assuming it vanishes at the boundary), giving \(\frac{d}{dt}\int\rho\,dx=0\). Applied to a probability density this means "the total probability is always 1" — the reason the normalization condition of the probability article is automatically maintained under the dynamics.

Relation to the transport equation. If \(v\) is constant and \(\nabla\cdot v=0\), then \(\nabla\cdot(\rho v)=v\cdot\nabla\rho\) and it reduces to the transport equation of the previous section. Conversely, if \(v\) varies with position, an extra term appears compressing or rarefying the density.

In machine learning: this is the bridge. When particles move according to the ODE \(\dot{x}=v(x)\), the density of those particles satisfies exactly this continuity equation. That is, the continuity equation is the dictionary joining "the language of trajectories" to "the language of densities." We confirm this directly in "From Particle to Density Dynamics" three sections from here, and the Probability Flow ODE of the final section uses the relationship in reverse — specifying a desired density change and finding the velocity field that produces it.

\[\begin{aligned} \frac{\partial \rho}{\partial t} + \nabla\cdot(\rho v) &= 0 \\ j &= \rho v \quad \text{(flux)} \\ \frac{d}{dt}\int \rho\,dx &= 0 \quad \text{(total mass conserved)} \end{aligned}\]

Play with different velocity fields. The density (blue) is carried along the velocity field, piling up where the flow converges and thinning where it diverges. The bar below is the total mass — check that it stays at 1 however much the shape changes.

Definition. \(\frac{\partial u}{\partial t}=\alpha\,\Delta u\). If \(u\) is temperature it is called the heat equation; if it is a concentration or a probability density, the diffusion equation — but mathematically they are exactly the same equation. Here \(\alpha>0\) is the diffusion coefficient.

What the Laplacian does. In the previous article the Laplacian was "how much lower this point's value is than the average around it." So this equation is a direct transcription of "rise if below your surroundings, fall if above them." Bumps are shaved away and dips filled in, and given enough time everything becomes flat.

The fundamental solution is a Gaussian. If the initial condition is a delta function concentrated at a point, the solution is \(u(t,x)=\frac{1}{\sqrt{4\pi\alpha t}}e^{-x^2/4\alpha t}\), a Gaussian whose variance grows linearly as \(2\alpha t\). The fact that the width spreads in proportion to \(\sqrt{t}\) reappears unchanged when we meet Brownian motion later — not a coincidence but two faces of the same phenomenon.

Why the fine ripples vanish first. Recalling that the Laplacian was "the difference from the surrounding average" makes the reason immediately visible. At the same amplitude, a tightly packed ripple travels from crest to trough over a shorter distance and so bends far more sharply, meaning a larger difference from the surrounding average. A larger difference is shaved away that much faster, so fine ripples are blurred out in an instant while only gentle large undulations survive. Confirming this intuition with formulas: the second derivative of a component \(\sin kx\) with wavenumber \(k\) is \(-k^2\sin kx\), so the decay rate is proportional to \(k^2\) and it dies as \(e^{-\alpha k^2 t}\) — double the frequency and the decay is four times as fast. Applying a Gaussian blur to a photograph, where fine texture goes first while large blocky forms remain, is the same phenomenon.

It is irreversible. Diffusion destroys information. The higher the frequency the faster it disappears (as \(e^{-\alpha k^2 t}\) in the wavenumber \(k\)), so after only a little time the fine structure is unrecoverable. The time-reversed equation \(u_t=-\alpha\Delta u\) blows up in the high frequencies and cannot be handled numerically. And yet diffusion models accomplish exactly this "reverse diffusion" — how that is possible is the story of the second half of this article (to give the answer away: because they learn the score, extra information about the data distribution).

In machine learning. The forward process of a diffusion model is blurring the data distribution according to this equation until it becomes a Gaussian. A Gaussian blur in image processing is the heat equation solved for a short time, and on a graph, diffusion by the graph Laplacian appears as message passing in a GNN and as the over-smoothing phenomenon.

\[\begin{aligned} \frac{\partial u}{\partial t} &= \alpha\,\Delta u \\ u(t,x) &= \tfrac{1}{\sqrt{4\pi\alpha t}}e^{-x^2/(4\alpha t)} \quad \text{(point initial condition)} \\ \hat{u}(t,k) &= \hat{u}_0(k)\,e^{-\alpha k^2 t} \quad \text{(high frequencies die first)} \end{aligned}\]

Play and the sharp peaks collapse first, followed by the broad structures. In the frequency graph below you can see the high frequencies dying exponentially faster — this is why diffusion is irreversible, and why rewinding does not restore the original shape.

Two descriptions. The Lagrangian view follows the particles: \(\dot{x}_i=v(x_i,t)\). The Eulerian view stands at a fixed location and counts what passes: \(\partial_t\rho+\nabla\cdot(\rho v)=0\). They describe the same physics, but their unknowns are trajectories and a density respectively.

The translation rule. Moving \(N\) particles along a velocity field \(v\) makes their empirical distribution converge, as \(N\to\infty\), to a solution of the continuity equation. Conversely, given a density, draw particles from it and let them flow. In the demo you can see the upper point cloud (particles) and the lower curve (density) moving together under the same velocity field and agreeing with each other.

Why both languages are needed. A density holds the complete information of the system but is impossible to store or compute in high dimensions (the same curse-of-dimensionality problem as a joint distribution). Particles can be handled even in high dimensions but are a finite sample and so carry noise. Building the theory in densities and computing with particles is the strategy in practice, and diffusion models work exactly that way — justified by Fokker–Planck and executed on sample paths.

Even deterministic motion spreads a density. Even if individual particles move without any randomness, the density spreads or gathers whenever the velocity field has divergence. That is, the diffusion of a density does not necessarily come from randomness. This observation is the seed of the final section's Probability Flow ODE — it is why the same density evolution as a random SDE can be reproduced by a deterministic ODE.

In machine learning. Generative models use this translation in both directions. Training deals with densities (likelihoods, scores), while sampling flows particles (a single path). Flow Matching first fixes a desired density path and then learns by regression the velocity field realizing it, and the correspondence of this section is exactly what makes that possible.

\[\begin{aligned} \dot{x}_i &= v(x_i, t) &&\text{(particles, Lagrangian)} \\ \frac{\partial\rho}{\partial t} + \nabla\cdot(\rho v) &= 0 &&\text{(density, Eulerian)} \\ \hat\rho_N = \tfrac{1}{N}\textstyle\sum_i \delta_{x_i} &\;\xrightarrow{N\to\infty}\; \rho \end{aligned}\]

Increase the number of particles N and play. The histogram of the upper particles gradually overlays the solution of the continuity equation below (orange). Lumpy when N is small and smoother as it grows — the convergence seen in the probability article's sampling section, meeting dynamics here.