Structured Simulators and Neural Operators

Suppose you are building a world model for a floating robot arm in a turbulent fluid. An unconstrained sequence model might learn to predict the next frame from video, but it may require substantial task data, can fail under distribution shift, and provides no architectural guarantee of momentum conservation or finite-speed propagation. Regularities such as angular-momentum balance and delayed pressure-wave propagation must then be inferred from observations. If you can instead encode rigid-body equations for the arm, Navier–Stokes equations for the fluid, and their coupling, valid structure can improve data efficiency, extrapolation, and interpretability. The tradeoff is that an authored model cannot represent an omitted effect by itself. Classical solvers may expose analytic, automatic, adjoint, or finite-difference sensitivities, but a solver that provides only forward outputs does not directly supply the derivatives needed for learning and planning.

This chapter is about the middle ground. We study four families of methods that inject physical structure into world models while keeping them learnable:

• Differentiable simulators: hand-built physics engines whose outputs can be backpropagated through, so the same machinery that trains neural networks can identify parameters and optimize plans. The physics is authored by hand; learning touches the parameters, the plans, and any neural components bolted on top.
• Graph-based physical models: learned dynamics whose architecture mirrors the fact that physical influence is local and relational, built from interactions between neighboring entities. The existence of an interaction law is assumed; its form is learned.
• Neural differential equations: neural networks that define a continuous-time vector field, so learned dynamics can be integrated, analyzed, and constrained with the tools of dynamical systems. The existence of smooth continuous-time flow is assumed; the flow itself is learned.
• Operator learning: neural architectures whose inputs and outputs are entire fields, learning maps from functions to functions for PDE-governed worlds like fluids, weather, and heat. The function-space character of the problem is assumed; the solution operator is learned.

These four families form a spectrum from "physics written by hand, learning only parameters" to "physics learned end to end, structure only in the architecture." Each step along that spectrum trades one kind of prior knowledge for another: the differentiable simulator assumes the equations, the graph model assumes locality and shared interaction laws, the neural ODE assumes a sufficiently regular continuous-time flow, and the operator network represents fields while still encoding architectural priors such as bases, kernels, geometry, and local or nonlocal coupling. Where you sit on that spectrum is one of the most consequential design decisions in a physical world model, and the right answer depends on how much data you have, how well you know the physics, and whether the model will be used to predict or to act. A model that only needs to forecast can tolerate structural compromises that a model used for planning cannot, because a planner can exploit inaccuracies that are reachable under its actions and favorable to its objective; not every inconsistency is reachable or useful.

Throughout Part I we argued that a world model is a learned surrogate for the environment's transition dynamics. Part III showed that the representation of state matters enormously: object-centric and relational state (Part III, Chapter 5) made latents more controllable and composable. This chapter makes the complementary point about the dynamics: the function that advances state through time can itself carry structure, and that structure is a form of prior knowledge you can spend instead of data. The two ideas compose rather than compete. A well-structured state makes the transition easier to learn; a well-structured transition makes the state representation's structure pay off, because compositional state is only useful if the dynamics respect that composition. An object-centric latent whose dynamics are a dense matrix over the flattened vector quietly destroys the decomposition the representation worked so hard to find.

The physical world has strong regularities that a generic function approximator must rediscover from scratch, and each of these regularities, once stated, suggests a family of architectures:

• Locality: many hyperbolic systems and explicitly local interactions propagate influence through neighboring regions at finite speed. Sound takes time to travel, and a disturbance at one end of a beam is not felt at the other until material waves carry it there. This motivates graph architectures with local message passing. The scope matters: diffusion has nonzero mathematical influence at arbitrary distance for positive time, and elliptic solution maps can depend globally on boundary data even when their differential operators are local.
• Relationality: many forces arise from interactions between pairs of entities, with the same law applied to each pair. Gravity applies the same pairwise law everywhere; its force depends on the bodies' masses and relative displacement, not arbitrary identities. This is what licenses sharing one learned interaction function across all entity pairs.
• Continuity and smoothness: states evolve along trajectories, not in arbitrary jumps, at least away from impacts and contacts. This is what licenses differential-equation models and integration as the mechanism of prediction.
• Conservation: in a closed, unforced system, conserved quantities such as mass, momentum, and energy restrict the space of plausible futures. Boundaries, forcing, dissipation, and exchanges with the environment introduce balance terms, so an apparent gain or loss must be interpreted with the system boundary made explicit.

When you hard-code these regularities, the model is reliable only where the encoded physics is valid. Learning without them offers more flexibility but may require more data and may generalize poorly. The methods in this chapter let you learn within a structured hypothesis space, so the data can focus on what you do not know. This is the same economy that shaped the identifiability discussion in Part II: valid structural assumptions shrink the hypothesis space and can reduce the observations needed for identification and improve generalization.

One caution before we begin, because it echoes a theme from Part I: a structured model that predicts well in open loop is not automatically a good model for control. Structure improves the trajectories a model produces, but planning exercises the model differently than prediction does: it pushes the model along action sequences it was never evaluated on, multiplies small errors over long horizons, and depends on gradients that can be wrong even when predictions look right. We will return to this distinction repeatedly, and it becomes central when these models drive planning in Part VII.

A simulator is a world model you wrote yourself. It takes a state $x_t$ and an action $a_t$ and produces a next state through hand-coded physics:

where:

• $x_t$, $a_t$: the state and the action at time $t$
• $x_{t+1}$: the next state produced by the hand-coded physics
• $f$: the hand-coded transition function that advances the state one step
• $\theta$: the physical parameters, collecting masses, lengths, friction coefficients, stiffnesses

The distinction between the structure of $f$ and the values of $\theta$ is the whole design. The structure encodes your physical understanding: which forces act, how they compose, what integration scheme advances the state. The parameters encode what you do not know precisely: the true mass of the payload, the actual coefficient of friction on the table, the real stiffness of a contact. Many legacy or externally wrapped simulators expose only forward outputs, even though classical sensitivity machinery may exist. A differentiable simulator additionally exposes derivatives of the implemented computation through analytic, automatic, implicit, adjoint, or custom derivative rules. For the rollout above, the quantities of interest are the following partial derivatives:

where:

• $T$: the rollout horizon, the number of simulation steps taken from the initial state
• $x_T$: the simulator state after $T$ steps, the output whose sensitivity we compute
• $\partial x_T / \partial \theta$: how the final state changes with each physical parameter
• $\partial x_T / \partial x_0$: how the final state changes with the initial condition
• $\partial x_T / \partial a_{0:T-1}$: how the final state changes with the entire sequence of actions applied over the rollout

Each of these three gradients powers a different capability. The parameter gradient drives system identification, the initial-condition gradient drives sensitivity analysis and state estimation, and the action gradient drives trajectory optimization. One derivative mechanism enables all three, which makes differentiability valuable for a world model.

Why does this matter? Three reasons dominate.

System identification becomes gradient descent. Part II, Chapter 3 treated system identification as an optimization problem and surveyed gradient-free and linear methods. Those methods are workhorses, but they can scale poorly with dimension and nonlinear dynamics. When the simulator is differentiable, identification of $\theta$ reduces to minimizing a loss $L(x_{1:T}^{\text{data}}; \theta)$ by backpropagation. The method avoids choosing parameter updates from forward evaluations alone; each iteration still simulates the trajectory, compares it with data, backpropagates through the rollout, and updates the parameters. The gradient combines information from each step of the rollout.

Planning becomes gradient descent. If you can differentiate the accumulated reward or a target cost through the rollout, you can optimize an action sequence directly. This is the modern descendant of shooting methods from classical optimal control (Part II, Chapter 4), and we will do exactly this in the worked example. Automatic differentiation replaces much of the system-specific derivative transcription with loss.backward(), provided every operation on the executed path has a valid derivative. The resulting gradient is exact for the implemented computation, not necessarily for omitted or idealized physics.

Learning can wrap around the simulator. Gradients flowing through simulation let you fit residual neural networks on top of physics, fit simulator parameters jointly with policy parameters, or backpropagate a learned perception module into a physical model. This pattern appears in some hybrid differentiable-simulation pipelines: the simulator supplies a structured backbone, while learned components model residual effects and end-to-end derivatives coordinate the trainable pieces. It is one design option for sim-to-real work, not a description of every sim-to-real pipeline.

The rollout is a deep computation graph: applying $f$ once is one layer, so a horizon-$T$ rollout is a $T$-layer network. Everything you know about training deep networks applies, with the twist that the layers are physics steps and the depth is dictated by the horizon you care about, not by an architecture choice. There are two standard ways to get gradients through it.

Unrolled backpropagation treats the whole trajectory as a feedforward graph and applies reverse-mode autodiff. This is what torch.autograd does automatically. It is exact and simple, and because reverse-mode autodiff computes the gradient of a scalar loss with respect to all inputs in a single backward pass whose cost is comparable to the forward pass, it is remarkably efficient per unit of information. Its weakness is memory: without checkpointing or recomputation, it stores intermediate states and memory grows linearly with $T$. For a long rollout of a high-dimensional state, memory can become the binding constraint, depending on state size, per-step compute, hardware, and the checkpointing or recomputation strategy.

The adjoint method avoids storing the trajectory by solving an augmented system backward in time. For continuous dynamics $\dot{x} = f(x, t; \theta)$ and a terminal cost $J = \Phi(x_T)$, which measures the terminal state, define the adjoint state $\lambda(t)$ by

where:

• $\lambda(t)$: the adjoint state, a costate variable that accumulates sensitivity information as it is integrated backward in time from $t = T$ to $t = 0$
• $\partial f / \partial x$: the Jacobian of the vector field with respect to the state, evaluated along the forward trajectory
• $\lambda(T)$: the terminal condition, equal to the gradient of the terminal cost $\Phi$ with respect to the final state $x_T$

and then the gradient is an integral along the backward pass:

where:

• $dJ/d\theta$: the gradient of the terminal cost with respect to the parameters
• $\lambda^\top$: the adjoint state (transposed), propagated backward from $T$ to $0$ along the reverse integration
• $\partial f / \partial \theta$: the Jacobian of the vector field with respect to the parameters
• The integral runs from $T$ down to $0$ because the adjoint system is solved backward in time; the minus sign comes from that reversed direction of integration

The intuition behind the adjoint state is worth pausing on: $\lambda(t)$ answers the question "if you perturb the state slightly at time $t$, how does the final cost change?" It starts at $T$ holding exactly the cost's sensitivity to the terminal state, and as it integrates backward, each step multiplies it by the transpose Jacobian of the dynamics, which converts sensitivity one step later into sensitivity one step earlier. The parameter gradient then follows by accumulating, along the backward sweep, how each parameter perturbation at each moment feeds into that running sensitivity. In the ideal continuous-adjoint formulation, storage can be independent of the number of accepted steps when the trajectory can be reconstructed. Practical implementations retain solver workspace or checkpoints and trade storage for recomputation. There is a pleasing unity here: discrete unrolled backpropagation is in fact an adjoint method for the discretized system; the distinction is practical, not fundamental. Which side wins depends on state dimension, horizon, numerical stability, and whether the integrator can reconstruct or restart the trajectory reliably.

Physics is not everywhere smooth, and this is the central engineering difficulty. The clean calculus of the previous section assumes $f$ is differentiable, and the most interesting parts of physical interaction are precisely where it is not.

• Contacts and impacts create hybrid maps whose event activation, timing, contact set, and normal can change non-smoothly with state or parameters. Within a fixed impact mode, a restitution map such as $v_n^+ = -e v_n^-$ is linear and continuous in the incoming normal velocity; the derivative difficulty comes from switching between modes and moving event times.
• Friction and stiction involve complementarity conditions (either the object sticks with force below a threshold, or it slips with force at the threshold), which are non-smooth by nature. The switch between sticking and slipping is a set-valued force law, and set-valued maps have no single Jacobian at the switching surface.
• Chaotic dynamics amplify perturbations along unstable directions approximately as $e^{\lambda T}$ over the regime where linearization applies. After several Lyapunov times, long-horizon prediction and gradient-based optimization can become badly conditioned. The derivative can still be mathematically correct; the problem is that it describes extreme sensitivity to initial conditions, parameters, and numerical error.

Practical responses depend on the engine and contact pipeline. Some differentiable simulators use compliant or softened contact, while others differentiate constraint, convex, or position-based formulations; frameworks such as DiffTaichi and Warp support more than one simulator formulation. Any surrogate changes the relationship between the implemented gradient and ideal rigid-contact physics. Against chaotic sensitivity, practitioners often differentiate over shorter windows and re-fit or re-plan frequently, exactly the receding-horizon pattern we will see in the worked example. An adjoint changes memory and computation, not the underlying sensitivity.

We can isolate one mechanism by computing parameter sensitivity in two linear systems. One transition contracts perturbations and the other amplifies them. This is not a chaos experiment; it shows how repeated Jacobian products can make an additive parameter sensitivity saturate or grow with rollout length.

The landscape is smooth and has a single clear sampled basin around the true value. Within $4 \le g \le 16$ at spacing 0.1, the scan shows no comparably low sampled alternative, but it does not certify a global optimum outside that interval or between samples. Missing physics can leave systematic residual error or distort the parameter landscape; whether this appears as a floor, a flat direction, or local minima depends on the data, loss, and compensating parameters. Inspecting the loss landscape and residual behavior is therefore useful evidence for identification, but not a proof that the model class is correct.

Differentiable simulators assume you know the equations. Graph-based physical models assume something much weaker: you know that the system is a network of local interactions. You do not know the force law, the material model, or even in detail what the entities are; you only know that influence is local, that it happens between pairs, and that the same rule applies everywhere. This is the dynamics-side twin of the object-centric representations from Part III, Chapter 5. There we argued that state should be decomposed into objects and relations; here we build the transition function out of the same decomposition. The payoff of matching the two is that the representation and the dynamics now agree about what the world is made of, and learning can happen in the natural coordinates of the problem.

The canonical particle example is the Graph Network Simulator (GNS) of Sanchez-Gonzalez and colleagues (2020). MeshGraphNets, introduced separately by Pfaff and colleagues, applies related graph machinery to mesh-based simulation. In these systems the "objects" are particles or mesh nodes, the edges connect spatial neighbors, and learned updates stand in for parts of a hand-designed simulator. The inductive bias has three parts:

• State is a set of nodes, each carrying local information (position, velocity history, material properties).
• Interactions are edges, and the dynamics of a node depend on its neighbors through learned message-passing functions, not on a global flattened vector.
• The update predicts a derivative (typically acceleration), which is integrated with a fixed physical scheme such as semi-implicit Euler.

The third part deserves emphasis because it imports a useful kinematic structure. Predicting acceleration and integrating it enforces a second-order update from acceleration to velocity to position. It does not by itself enforce Newton's second law, momentum conservation, or force symmetry; those require additional parameterization or constraints. A model that predicts positions directly has no comparable integration scaffold.

A message-passing round looks like

where:

• $h_i$: the feature vector (for example, position and velocity history) of node $i$ before the update
• $h_j$: the feature vector of neighboring node $j$
• $e_{ij}$: the edge features carrying relative geometry between the pair, for example $x_j - x_i$
• $m_{ij}$: the message sent from node $j$ to node $i$, computed by the learned interaction network $\phi$
• $\sum_j m_{ij}$: the sum of all incoming messages to node $i$, aggregating its neighborhood
• $\psi$: the learned update network that combines a node's own features with its aggregated messages
• $h_i'$: the updated node features after one round of message passing

Here $\phi$ and $\psi$ are small neural networks, $e_{ij}$ carries relative geometry (for example $x_j - x_i$), and the sum aggregates all incoming messages. Two design details carry most of the physical weight. First, if all positional inputs enter through relative displacement rather than absolute position, translating the whole system leaves the interaction input unchanged. Other features, such as masses, material types, or relative velocities, may still affect the interaction. Second, summation matches additive force composition when messages are parameterized as pairwise force contributions; it does not alone enforce antisymmetry or conservation. Stacking $M$ rounds lets information travel $M$ hops per step for a fixed graph.

The payoff of this bias shows up in two places. First, sample efficiency: the model must learn one interaction law rather than a separate rule for every configuration of the system. A pair of neighboring particles observed in one corner of the domain constrains the same function that governs particles elsewhere, so every pair observation in the training data informs the single shared law. Second, size and shape transfer: a model trained on a chain of 8 masses can drive a chain of 12, because nothing in its parameters is tied to the number 8. The learned function $\phi$ takes a pair of features and returns a message; it neither knows nor cares how many such pairs exist. A fixed-width dense model over a flattened state does not natively accept another entity count; supporting size variation requires padding or masking, a maximum-size parameterization, shared structure, or architectural redesign.

Let us demonstrate the core mechanism in its simplest possible form. We will simulate a chain of masses connected by springs with an unknown rest length, assume we know only that (a) forces are pairwise between neighbors, (b) the same law $\varphi$ applies to every pair, and (c) damping is linear and known, and then recover $\varphi$ from observed motion by linear regression. This strips the message-passing machinery down to its identifiable core: what survives when you remove the neural networks is the hypothesis "one shared pairwise law," and that hypothesis alone is enough to make the problem solvable from modest data.

The learned law rolls out with tiny error on both tested chain sizes. Because the same local law applies at both sizes, data from 8-mass chains transfers to this matched 12-mass case. Transfer to other sizes still depends on coverage of the relevant local configurations and dynamics. The error curves also illustrate compounding rollout error: small per-step mismatches accumulate over the horizon. Full GNS models retain shared local updates while replacing this linear regression with neural message passing. Their learned interactions remain trustworthy only where training data or structural constraints support them.

Both previous approaches keep an explicit integrator and learn either parameters or a force law. Neural differential equations go one step further back: the learned object is the vector field. There is no hand-written equation and no assumed interaction structure; there is only the commitment that state evolves as a continuous-time flow, and a neural network is installed as the right-hand side of the flow equation. A neural ODE (Chen and colleagues, 2018) defines state evolution by

where:

• $z(t)$: the state of the system at continuous time $t$
• $f_\psi$: a neural network, with parameters $\psi$, that defines the time derivative of the state

Trajectories are produced by a compatible ODE solver, and the model's job is to learn the right-hand side from data. The change of perspective is subtle but consequential: in a differentiable simulator the integrator is an implementation detail of your physics, while here the integrator is a generic engine and the entire modeling content lives in the vector field. The commitment is that state evolution can be represented by a sufficiently regular continuous-time vector field; the regularity needed for existence, uniqueness, and differentiability depends on the model and solver.

The connection to the adjoint method is especially clean. The continuous-adjoint formulation can avoid storing every accepted solver step by reconstructing or recomputing the trajectory during the backward solve. Its idealized memory cost is independent of the number of stored steps, excluding solver workspace and checkpoints, but numerical error and backward instability can make its gradients differ from differentiating the actual discrete solver. The appropriate choice is therefore implementation- and problem-dependent.

What does continuous time buy you, concretely?

• Decoupling the model from the integrator. Because the learned object is a vector field rather than a discrete map, it can be evaluated with another compatible solver or time grid. Accuracy and stability can still change with the solver, tolerance, and sampling regime, so deployment at a new resolution requires validation rather than an assumption of invariance.
• Irregular sampling. Observation times need not be uniform. The solver simply integrates between whatever times data arrives, which is awkward for fixed-window discrete models and natural here. This connects directly to the temporal state and memory questions of Part III, Chapter 4. Real sensors do not tick at metronomic intervals, and a model class that treats the timestamps as part of the problem rather than a nuisance is materially better suited to them.
• Dynamical-systems tools. You can compute fixed points, Jacobians, Lyapunov exponents, or conservation properties of $f_\psi$, and you can constrain it: divergence-free fields for incompressible flow, Hamiltonian structure for conservative systems (the Hamiltonian and Lagrangian neural network line of work). This is a two-way street: the tools analyze what was learned, and the constraints shape what can be learned. An exact autonomous Hamiltonian vector field conserves its Hamiltonian along the continuous flow, but a numerical rollout can still drift unless the integrator and step-size regime preserve the relevant structure accurately.

The main training design choice is one-step supervision versus unrolled supervision, and it is the same open-loop-versus-compounding issue we have flagged before. Supervising $f_\psi$ to match $\frac{z_{t+1} - z_t}{\Delta t}$ at every pair of adjacent states is easy and stable, but it is exactly the one-step trap from our discussion of rollout error: tiny one-step biases integrate into large trajectory biases. A vector field that is locally correct to within a percent can still point systematically off-trajectory, and over a hundred integration steps those small angular errors push the solution onto a path the true system never visits. Training through short unrolled segments pushes the model toward trajectories that stay on track, and is usually worth the extra cost.

Let us learn a Van der Pol oscillator from a single trajectory. Its dynamics,

where:

• $z_1, z_2$: the two state variables, position-like and velocity-like respectively
• $\mu$: the nonlinearity parameter; with $\mu > 0$ the system exhibits a self-sustained limit cycle, and here $\mu = 1.5$
• $\mu(1 - z_1^2) z_2$: the state-dependent damping term, which removes energy for large $|z_1|$ and injects energy near the origin, producing the limit cycle

produce a self-sustained limit cycle: trajectories spiral out from near the origin and settle onto a closed loop. The Van der Pol system is a demanding target for exactly the reason it is a good test: its long-term behavior is a geometric object, a closed attracting loop, and reproducing it requires getting the energy balance right everywhere around the cycle, not just on average. A learned vector field that injects slightly too much energy will spiral outward forever; one that removes slightly too much will collapse to the origin. The learned model sees only states, never the equations.

The learned rollout traces a loop close to the visited limit cycle. From segments of one clean simulated trajectory, this run learns a vector field whose long rollout visually approximates the oscillation. The figure does not by itself establish global vector-field recovery or precise amplitude and period errors. Two failure modes are worth knowing about before you rely on this pattern in earnest.

Stiffness. If the true dynamics have fast transients relative to the slow behavior of interest, explicit solvers need tiny steps, and training through them is expensive. The problem is not merely slow training: an adaptive solver may take very different numbers of internal steps for nearby inputs, which makes the training landscape itself jagged. Recognizing stiffness and switching to implicit solvers is a real skill, and the adjoint method helps because it does not need to store the fine-grained forward trajectory.

Drift and non-conservation. Nothing in the architecture prevents the learned vector field from slowly spiraling inward or outward where the true system conserves energy. The Van der Pol example survives this because the limit cycle is attracting, so small energy errors are corrected by the dynamics themselves; a conservative system with no attractor gives the learned model nowhere to hide, and unbounded drift accumulates. Structural variants address this: Hamiltonian and Lagrangian neural networks build conservation into the parameterization, and augmented neural ODEs add extra dimensions so the model can express trajectories that a fixed-dimensional field cannot fit without distortion. Uncertainty about where the model is trustworthy, rather than a single confident trajectory, is the deeper fix, and that is the subject of the next chapter.

Neural ODEs are also a conceptual bridge to Part V. Residual recurrent state updates can sometimes be interpreted as discretizations of continuous-time latent dynamics, with an encoder mapping observations into that state. Recurrent state-space models in general are broader: they may be gated, stochastic, discontinuous, noninvertible, or otherwise not embeddable in a same-dimensional ODE flow.

Everything so far concerns systems represented with finite-dimensional state. Continuum fields, such as fluid velocity, material temperature, or ocean pressure, are functions governed by partial differential equations. A numerical discretization with 10 million unknowns is a large finite-dimensional ODE or DAE approximation of that continuum problem. A fixed-grid architecture may couple its learned representation to one resolution or geometry, while convolutional, graph, interpolation, multigrid, and operator constructions can share parameters across compatible discretizations. Whether they transfer depends on coordinates, normalization, geometry handling, numerical resolution, and the training distribution.

Operator learning reframes the learning target. For a PDE with solution $u(x, t)$, the object of interest is often not one trajectory but the solution operator: the map

which takes an entire initial condition to the entire solution at time $T$. A neural operator is designed to approximate this function-space map from sampled representations. Appropriately constructed operators can share parameters across discretizations and sometimes evaluate zero-shot at new resolutions. That capability is not automatic: standard FFT-based FNOs assume regular grids, while new geometries need compatible coordinates, meshes, or geometry-aware variants and remain distribution-dependent.

Two foundational architectures illustrate different answers to the question "how does a network eat a function?"

DeepONet (Lu and colleagues, 2019) factorizes the output function evaluation as a learned inner product:

where:

• $G(u)(y)$: the operator output at the query location $y$ for the input function $u$
• $b_k(u)$: the $k$-th coefficient the branch network produces from the entire input function $u$
• $t_k(y)$: the $k$-th feature the trunk network produces at the query location $y$
• $p$: the number of coefficient-feature pairs, which sets the expressiveness of the factorization

A branch network encodes sampled values of the input function $u$ into coefficients $b_k$, and a trunk network encodes each supported query location $y$ into features $t_k$. Their dot product lets one trained model evaluate new query points within the represented domain, subject to the sensor layout and training distribution. The factorization resembles a learned basis expansion: the trunk provides basis-like spatial functions and the branch provides input-dependent coefficients. Universal approximation results apply under compactness and continuity assumptions; they are existence results, not finite-data generalization guarantees.

The Fourier Neural Operator (FNO, Li and colleagues, 2020) builds the nonlocality into a single layer. Its core is the spectral convolution:

where:

• $v$: the input field (for example, the feature values on the grid at the current layer)
• $\mathcal{F}$: the Fourier transform, mapping the field to its frequency modes
• $R$: a learnable complex weight applied to each frequency mode, with the highest modes truncated
• $\mathcal{F}^{-1}$: the inverse Fourier transform, mapping the weighted modes back to the grid
• $(K v)(x)$: the output of the spectral convolution at grid location $x$

Here $\mathcal{F}$ is the Fourier transform, $R$ is a learnable complex weight applied to each frequency mode (with the highest modes truncated), and $\mathcal{F}^{-1}$ transforms back. Multiplication by a learned $R$ in frequency space is a convolution on the grid, and unlike a local CNN kernel it can couple every point to every other point in one operation. This global path is useful for solution operators with nonlocal dependencies, including positive-time diffusion on the periodic domain used below; not every PDE or forecast horizon requires all-to-all influence. A full FNO stacks several layers of the form

where:

• $v_l(x)$: the layer-$l$ feature value at grid location $x$
• $W$: a local linear map applied pointwise at each location
• $(K v_l)(x)$: the global spectral path from the previous equation, coupling every location to every other
• $\sigma$: a pointwise nonlinearity applied after summing the local and global paths

The layer combines the global spectral path with a local linear path and a pointwise nonlinearity. Nonlinear PDEs couple modes, while a purely linear spectral map remains confined to linear transformations regardless of depth. The pointwise nonlinearity lets successive spectral layers compose nonlinear operators.

Why is the spectral view so natural for PDEs? Because many linear PDE solution operators are diagonal in the Fourier basis. The heat equation $u_t = \nu u_{xx}$ on a periodic domain multiplies each wavenumber-$k$ Fourier coefficient by $\exp(-\nu k^2 T)$. A spectral architecture can represent this multiplier exactly on its retained modes; an architecture that looks at one grid point at a time cannot express the full map for arbitrary fields, because $u(x,T)$ depends on the entire initial profile. When you choose an architecture for a field-valued world, you decide which dependency structures are available to the model.

Weather forecasting provides large-scale examples. In their 2022--2023 evaluations, FourCastNet, based on an adaptive Fourier operator, matched then-current IFS performance on selected short-lead variables, while GraphCast, a multimesh graph model, beat then-current HRES on most specified verification targets. These were benchmark-, lead-time-, hardware-, and baseline-specific results rather than universal dominance. FourCastNet reported a week-long forecast in under two seconds on its stated accelerator setup; GraphCast reported a ten-day forecast in under 60 seconds on a TPU. Both timings describe inference and exclude training and operational pipeline costs.

Let us see the argument play out on the heat equation. We generate discrete spectral solutions, then train two models on the same data: a minimal spectral operator and a pointwise MLP that sees only $u_0(x)$ at each location. The spectral architecture can represent the heat multiplier on its retained modes; modes above its cutoff are present but strongly attenuated in the targets. Both models use the same data and number of training steps. The spectral model has 2,073 trainable real parameters and the pointwise MLP has 4,353, so this comparison isolates the dependency structure more clearly than it isolates parameter count.

The contrast is stark. The spectral operator tracks the target closely, while the larger pointwise MLP does not. The failure is an expressiveness problem imposed by the dependency structure: no function of the single number $u_0(x)$ can equal a value that depends on the entire field for arbitrary initial conditions. Diffusion damps high wavenumbers, so by time $T$ the output at one point is a weighted average of values across the periodic domain. A pointwise map has no mechanism for that averaging. This fixed-grid experiment isolates nonlocal representability; it does not test transfer to another resolution or geometry.

Everything above was about prediction quality. This worked example turns a structured model toward decision-making, which is where world models earn their keep, and it demonstrates the difference between open-loop prediction and closed-loop usefulness that runs through this book. The example is intentionally simple, because its purpose is not to showcase a hard control problem but to expose, in a setting where the answer is known, exactly what differentiability contributes to planning and exactly where the resulting plans are fragile.

The task: a damped pendulum starts at angle $\theta_0 = 1.2$ and we want the final angle near $\theta = 0$ after 2.4 seconds, using a bounded torque $|u| \le 2$. We use normalized units with effective moment of inertia $I=1$, so torque enters the angular acceleration as $+u$; otherwise the dynamics contain $+u/I$. We know the physics exactly, so the differentiable simulator is the world model, and the plan is simply a sequence of torques found by gradient descent through the rollout:

where:

• $u_{0:T-1}$: the sequence of $T$ torque actions being optimized, the decision variables of the problem
• $\theta_T$: the pendulum angle after $T$ steps under the chosen action sequence, which we want near zero
• $\rho$: the effort weight penalizing control energy, set to $0.01$ in the experiment below
• $\frac{1}{T}\sum_t u_t^2$: the mean squared normalized torque over the horizon, discouraging aggressive actuation
• $x_{t+1} = f(x_t, u_t)$: the differentiable simulator dynamics as a constraint on every step
• $u_{\max}$: the torque bound, with $|u_t| \le u_{\max}$ enforced via a tanh squashing in the implementation

This is a shooting method from Part II, Chapter 4, with autodiff replacing hand-derived adjoints. The constraint is handled implicitly: rather than enforcing the dynamics as constraints in a constrained optimizer, we simulate them forward and optimize only over the action sequence, which is what "shooting" means. The dynamics are satisfied by construction inside the simulator; the price is that the parameterization of the decision variables, a flat sequence of torques, gives the optimizer no mechanism to correct early mistakes except through the actions themselves.

Gradient descent through the simulator finds a torque sequence that reaches the target angle at the terminal time. The pendulum is still moving, because the objective penalizes terminal angle but not terminal angular velocity or dwell time. Notice what made the optimization work: differentiability of the implemented physics supplied action sensitivities, differentiability through time composed them across 120 steps, and differentiable squashing enforced the torque bound at every iterate. These gradients are correct for the implemented model; mismatch between that model and reality can still mislead the planner.

This plan is open loop: we optimized one torque sequence against the model and executed it. If the real pendulum's friction differed from the model's, or a nudge arrived mid-swing, the plan would drift and the realized terminal error would differ from the prediction. The plan is most fragile where the optimized trajectory passes through sensitive regions of state space, yet a single deterministic rollout does not identify those regions. The standard remedies are to re-plan from the measured state every few steps (receding-horizon control, developed in Part VII) and to model residual uncertainty explicitly. Re-planning replaces one long plan with a chain of shorter optimizations at the cost of solving a problem at each control step. A model accurate enough to optimize through is not automatically safe to trust for the full horizon. This open-loop-versus-closed-loop gap is a central evaluation criterion for the architectures in Part V.

The four families form a design spectrum, and the choice among them is driven by what you know, what you have, and what the model is for. Reading the table by columns is instructive: the "hand-designed" column generally shrinks as you move down, while the "sweet spot" column tracks how much physics you can commit to in advance. The data-hunger labels are qualitative heuristics, not controlled cross-family measurements; they depend on the task, target error, architecture, data quality, and evaluation regime.

A few rules of thumb distilled from the chapter:

• Differentiate what you know, learn what you do not. If the equations are right up to parameters, use a differentiable simulator and fit the parameters. If only the interaction structure is known, learn inside a graph model. Reserve learning capacity for unknown components. Every valid equation you supply is one less relation the data must identify.
• Match the architecture's bias to the physics' structure. Pairwise locality calls for graphs, temporal smoothness calls for differential equations, and spatial nonlocality calls for operators. A mismatched bias can waste capacity or exclude the target dependency entirely, as the pointwise heat baseline demonstrates.
• Prefer hypothesis classes that contain the target dependency structure and can approximate the dynamics to the fidelity the task requires. The pointwise MLP failed the heat equation not for lack of width but because it had no pathway for cross-location dependence. Extra width cannot repair a missing dependency pathway.
• Hybrid physics-plus-residual models are useful when core physics is trustworthy and systematic residual effects are learnable. Physics-informed neural networks likewise add differential-equation residuals and boundary or initial-condition information to the objective, although some variants enforce conditions directly. A persistent, validated residual can flag where a mechanistic model disagrees with data, but attributing it to missing physics requires ruling out observation noise, state-estimation bias, discretization error, parameter error, distribution shift, and failed optimization.

The deepest limitation of every method in this chapter is that structure is a wager on where the model will be used. A differentiable simulator is accurate only inside the regime its equations describe; a graph network's learned interaction law assumes coverage of relevant pair configurations; an operator trained on one family of initial conditions may degrade outside it. Agents that optimize through these models can find the seams. A policy may exploit compliant contact in ways a real robot cannot, while an operator-based planner may seek inputs where prediction error flatters its objective. Structure narrows the space of plausible failures, but it does not eliminate model exploitation. Part XII returns to these reliability questions.

The second systemic limitation concerns evaluation. These families are often evaluated on held-out predictive error---trajectory or rollout error for dynamical models, and field or function error for operator models---but predictive error alone does not establish usefulness for control. As we saw in the worked example, open-loop accuracy and closed-loop usefulness can come apart. A simulator whose parameters are slightly wrong may produce visually indistinguishable trajectories while inverting the sign of a gradient that planning depends on; a graph model with excellent one-step accuracy may drift over horizons that matter for a task. Structured models earn their complexity budget only when the structure they encode is the structure the decision problem exercises, and evaluating that requires task-level protocols, not just predictive-error metrics. In practice this means evaluating plans executed under model mismatch, not just predictions compared to data, and it means reporting where in state space the model is trustworthy, not just how wrong it is on average.

The impact of these methods has been substantial. Differentiable simulators connect system identification, optimal control, and deep learning. Graph simulators have shown that learned local updates can model particle and mesh systems, and GraphCast-style models scale relational structure to planetary forecasting. Neural differential equations connect learned dynamics to dynamical-systems analysis, while operator learning targets reusable PDE solution maps. Each capability remains conditional on the encoded structure, numerical method, data regime, and deployment distribution.

One more frontier remains deliberately open in this chapter. Every model here produced a deterministic prediction. Physical prediction is uncertain because chaotic dynamics amplify initial-state uncertainty, contact conditions and parameters are imperfectly known or variable, and sensors are noisy; only some of these mechanisms are intrinsically stochastic. Representing aleatoric, epistemic, and observation uncertainty changes both training objectives and how a planner should consume predictions. A predictive distribution can support risk-aware decisions, provided it is calibrated in the relevant state-action region. That is the subject of the next chapter on stochasticity and uncertainty.

• Structured simulators inject physical knowledge at different points: equations (differentiable simulators), interaction topology (graph models), temporal structure (neural ODEs), or function-space structure (operators). Choosing the injection point is a core world-model design decision, because it fixes what the data must teach and what the architecture guarantees for free.
• Differentiable simulators turn identification and planning into gradient-based optimization, with unrolled backpropagation trading memory for simplicity and adjoints trading storage for recomputation or backward integration. Contact formulations are engine-specific, and chaotic sensitivity often motivates shorter differentiation and replanning horizons.
• Graph-based physical models encode locality and shared pairwise interaction laws. Parameter sharing can improve sample efficiency, while the spring-chain experiment demonstrates transfer from the matched 8-mass training setup to a 12-mass test. Interior motion only constrains differences of pairwise forces, so boundary information or priors are needed to pin down an absolute law. Learned interactions remain limited by their training distribution.
• Neural ODEs learn the vector field itself, decoupling the model from the integrator and enabling dynamical-systems analysis. Unrolled training mitigates one-step-to-trajectory drift; conservation requires structural commitments, and stiff systems demand solver sophistication.
• Operator learning targets maps between function spaces. Spectral architectures can represent nonlocal operators that pointwise models cannot. Compatible neural-operator designs may share parameters across discretizations, but the fixed-grid heat experiment tests nonlocal representability, not cross-grid transfer.
• Structure trades data for assumptions. Valid priors can reduce data needs and improve extrapolation within their regime, while incorrect priors can fail sharply. Hybrid residuals are useful corrections and potential diagnostics, but a large residual is not uniquely attributable to missing physics.
• Prediction quality is not decision quality. Gradients through a model can optimize plans that the same model's open-loop errors will break; receding-horizon re-planning and explicit uncertainty modeling are the standard responses, and uncertainty is where the next chapter begins.

Ready to test your understanding? Take this quick quiz to reinforce what you've learned about structured simulators and neural operators.