Dense Output

When solving ordinary differential equations (ODEs) numerically, we often get a series of solution points at discrete, unevenly spaced moments in time. An adaptive solver might take large steps when the solution is smooth and tiny steps when it changes rapidly. This raises a critical question: what is the solution's true behavior in the gaps between these points? Simply connecting the dots with straight lines is a naive approach that can introduce significant errors, miss crucial events, and fundamentally misrepresent the system's dynamics. This gap in our knowledge poses a major challenge for many scientific and engineering problems that require a continuous understanding of a system's evolution.

This article explores the elegant solution to this problem: ​​dense output​​. You will learn how modern ODE solvers, far from being simple point-to-point jumpers, can be designed to provide a high-fidelity, continuous representation of the solution at minimal extra cost. The following chapters will guide you through this powerful numerical method.

• ​​Principles and Mechanisms​​ will explain what dense output is, how it leverages the internal workings of a solver to generate a trustworthy continuous function, and why this approach is vastly superior to basic interpolation.

• ​​Applications and Interdisciplinary Connections​​ will demonstrate how dense output becomes an indispensable tool, enabling precise event detection, the rigorous analysis of chaotic systems, the solution of equations with memory, and even bridging the gap between continuous models and digital control systems.

Imagine you are trying to map the path of a satellite orbiting the Earth. You can't track it continuously; instead, you get a series of snapshots in time: its position here at 1:00 PM, its new position there at 1:03 PM, another over there at 1:04 PM, and so on. Your job is to describe its full, continuous trajectory. But the snapshots, taken by your numerical solver, are like stepping stones across a river—the time gaps between them are not equal. The solver, being clever, takes large steps when the satellite is coasting smoothly through empty space and tiny, careful steps when it's maneuvering or passing near another body.

This leaves us with a fundamental riddle: what happened in the gaps between the stones?

The most straightforward idea is to just connect the dots. If we have the satellite's position at time $t_n$ and the next at $t_{n+1}$, we can draw a straight line between them. This is called ​​linear interpolation​​. For many simple purposes, it's not a terrible guess. If the time steps are very, very small, the true path probably doesn't curve much, and a straight line is a decent stand-in.

But what if the steps aren't so small? What if during that interval, the satellite fired a thruster or swung by the Moon? The true path would have a significant curve, a bend that our straight line completely misses. A plot made of straight-line segments will look jagged and artificial, and more importantly, it can be misleading. If we wanted to know the exact moment the satellite crossed the Earth's equator, a straight-line approximation between a point in the northern hemisphere and a point in the southern hemisphere could give us a time that is quite wrong.

This is a deep problem in numerical mathematics. Simply forcing a curve through a set of points doesn't guarantee you've captured the true nature of the function between those points. In some cases, using more and more points to build a higher-degree polynomial can make the approximation worse, leading to wild oscillations that weren't there in the first place. The lesson is that naive interpolation is a dangerous game. We need a more principled approach.

This is where the magic of modern ODE solvers comes in. To decide how big a step to take from $t_n$ to $t_{n+1}$, the solver has to do a great deal of work. It doesn't just look at the beginning and end; it "probes" the dynamics inside the interval. For instance, a famous method like the Dormand-Prince 5(4) pair—the engine behind ode45 in many software packages—calculates the system's rate of change at several intermediate points within the step. These are called ​​stage derivatives​​.

For a long time, this intermediate information was treated as temporary scrap, used only to estimate the error and decide the next step size. Once the step was accepted, all that rich detail about the "shape" of the solution within the interval was thrown away.

The concept of ​​dense output​​ is born from a simple, brilliant insight: Don't throw that information away!

Instead of just keeping the endpoints $(t_n, \mathbf{y}_n)$ and $(t_{n+1}, \mathbf{y}_{n+1})$, a solver with dense output uses the "ghost" information—the internal stage derivatives it already calculated—to construct a high-quality interpolating polynomial, let's call it $\mathbf{y}_{interp}(t)$, that is valid across the entire interval $[t_n, t_{n+1}]$. This isn't a simple straight line. It's a smooth, continuous curve, typically a polynomial of degree 3, 4, or 5, that not only connects the two endpoints but also respects the curvature and flow of the solution inside the interval. This is possible because the design of the solver itself was optimized for this purpose. The celebrated Dormand-Prince method, for example, was developed not just for its accuracy and efficiency (achieved through clever tricks like FSAL, or "First Same As Last"), but also specifically to provide a high-quality continuous extension.

So, what makes this special? The crucial advantage of dense output is that the error of this interpolating polynomial is consistent with the error of the solver itself.

Let's be clear about what this means. When you set an error tolerance for an adaptive solver, you are essentially telling it, "I want the error at each stepping stone to be no more than this amount." The solver then adjusts its step sizes to meet that promise. But linear interpolation between those stones offers no such guarantee. The error between the stones could be much, much larger than your tolerance.

Dense output, however, extends the promise. Because the interpolant $\mathbf{y}_{interp}(t)$ is constructed using the same underlying mathematics that controls the step-wise error, it maintains the same order of accuracy across the entire interval. If your solver is a 5th-order method delivering a solution with a local error of, say, $10^{-6}$, the dense output guarantees that the value you get from $\mathbf{y}_{interp}(t)$ at any time $t$ between $t_n$ and $t_{n+1}$ also has an error of roughly that same order. You get a continuous function that you can trust to the same degree as the solver's discrete points.

This high-fidelity continuous solution is not just an academic nicety; it is a game-changer for a vast range of practical problems.

• ​​Finding Hidden Events:​​ Let's go back to the engineer from our problem who needed to find the exact times when a component of her oscillating system crossed zero. These events almost never happen to fall exactly on one of the solver's discrete time steps. With dense output, she has a high-accuracy function, $\mathbf{y}_{interp}(t)$. Finding the event is now a matter of finding the roots of this function (e.g., solving $y_{interp,1}(t^*) = 0$), a standard and accurate numerical task. Trying to do this with crude linear interpolation would yield much less precise results.

• ​​Smooth Visualizations:​​ If you want to create a smooth animation or a publication-quality plot, dense output is the right tool. You can query the interpolant at as many points as you need to generate a visually perfect curve, without the computational cost of forcing the solver to take tiny steps all the way.

• ​​Delay-Differential Equations (DDEs):​​ Some of the most interesting systems in biology, economics, and control theory have memory. The rate of change of the system now depends on its state at some time in the past. A simple model of a population $P(t)$ with a maturation delay $\tau$ might look like $\frac{dP}{dt} = r P(t) \left(1 - \frac{P(t-\tau)}{K}\right)$. To calculate the derivative at time $t$, the solver needs to know the value of $P(t-\tau)$. The point $t-\tau$ is almost certainly not a point where the solver has a stored "stepping stone". Dense output solves this problem elegantly. It provides a history of the solution that can be queried at any past time, making the numerical integration of DDEs possible.

In essence, dense output transforms an ODE solver from a machine that produces a discrete sequence of points into one that produces a true, functional representation of the solution. It is a beautiful example of how paying attention to the "scraps" of a computation and designing algorithms with a holistic view can unlock powerful new capabilities, turning a simple dot-to-dot plotter into a master artist.

Now that we have taken the engine of our ordinary differential equation (ODE) solver apart and seen how the gears mesh, it's time to take it for a drive. The real fun of a tool isn't just knowing how it works, but seeing the vast and unexpected landscapes it allows us to explore. A standard adaptive solver that only reports the solution at its chosen time steps is like teleporting through a beautiful countryside—you see the destinations, but you miss the entire journey. Our "dense output" is the magnificent bay window on our vehicle, allowing us to see everything that happens between the steps.

In the previous chapter, we learned that dense output provides a high-quality continuous function that approximates the solution across the intervals chosen by the solver. It doesn't just connect the dots with straight lines; it uses the dynamics of the system—the derivatives—to draw a smooth, accurate curve. This seemingly simple feature transforms a discrete-point generator into a powerful function approximator, and in doing so, it unlocks solutions to problems across a remarkable range of scientific disciplines.

Perhaps the most intuitive application of dense output is in finding the exact moment something interesting happens. Imagine you are simulating a rocket launch. You don't just want to know its position at $t=1, 2, 3$ seconds. You want to know the exact time it reaches its apogee (when its vertical velocity becomes zero), the precise moment of main engine cutoff, or the instant it crosses the Kármán line and officially enters space. These are "events," and they rarely occur exactly at the discrete time points your solver happens to choose.

If you only have the discrete solution points, you might see that the rocket's altitude was below 100 km at $t_k$ and above 100 km at $t_{k+1}$. You could try to guess the crossing time by drawing a straight line between the two points, but this is just an approximation. The rocket's true path is a curve, not a line, and this linear interpolation could be off by a significant amount. Worse, if the event happens and reverses itself entirely between two large steps—like a brief voltage spike—your solver might miss it completely!

This is where dense output, combined with a numerical root-finding algorithm, provides a spectacular solution. We define an "event function," $g(t)$, whose roots correspond to our event. For the Kármán line example, this would be $g(t) = \text{altitude}(t) - 100 \text{ km}$. The solver integrates as usual, but after each step from $t_k$ to $t_{k+1}$, it checks if the sign of $g(t)$ has changed. If it has, a root must lie within the interval. The solver then hands its dense output—the high-quality polynomial representing the solution in that interval—to a root-finder. This acts like a computational microscope, allowing us to zoom in on the interval and find the time $t^*$ where $g(t^*) = 0$ with a precision that matches the solver's own accuracy.

This capability is indispensable in control engineering. For instance, when analyzing the response of a system to a step input, a key metric is the "settling time"—the time it takes for the output to enter and remain within a certain percentage (say, 5%) of its final value. For a general, complex system that can't be solved with pen and paper, locating this precise moment is a classic event-finding problem, elegantly solved with dense output.

The need for precision becomes even more critical when we venture into the bewildering world of nonlinear dynamics and chaos. Here, tiny errors don't just lead to small inaccuracies; they can lead to conclusions that are qualitatively and fundamentally wrong. A key tool for understanding chaotic systems is the ​​Poincaré section​​, which acts like a stroboscope, freezing the motion of a system each time it passes through a specific plane in its state space. By plotting these intersection points, we can transform a tangled, impenetrable trajectory into a structured, often beautiful, pattern—the strange attractor.

But how do you record these intersection points? A naive approach might be to just take whichever computed point lies closest to the plane. Another approach might be to find two adjacent points on opposite sides of the plane and use linear interpolation to estimate the crossing point. Both, it turns out, are terrible ideas when using an adaptive step-size integrator.

The pitfall is subtle and profound. An adaptive solver is not a uniform observer. It takes small steps where the dynamics are complex and the trajectory is curving sharply, and it takes large, leisurely steps where the motion is smooth and predictable. Imagine trying to study the demographics of a bustling city by sending out a photographer who is instructed to take a photo every time they "feel something is changing." They would take dozens of photos of a thrilling street performance but might walk for an hour through quiet residential streets without taking a single shot. The resulting album would suggest the city is nothing but street performers!

This is exactly what happens with naive sampling of a chaotic system. By oversampling the "interesting" (high-curvature) regions and undersampling the "boring" (low-curvature) ones, you create a systematically biased picture of the attractor. The beautiful fractal structure you hope to see becomes warped and distorted, and any statistical properties you try to calculate will be wrong.

Dense output provides the only truly rigorous solution. It allows you to be a perfect, unbiased observer.

• For a ​​geometric section​​, where you want to find where the trajectory crosses the surface $h(\mathbf{x}) = 0$, you use the continuous output and a root-finder to find the exact intersection point for each crossing, just as in our event-finding example.
• For a ​​stroboscopic section​​ in a periodically forced system, where you need to sample at exact time intervals $t = T, 2T, 3T, \dots$, you use the dense output to evaluate the solution at these precise moments, regardless of where the solver's natural steps fell.

In doing so, you remove the systematic bias introduced by the adaptive stepping and recover the true geometry of the attractor in all its intricate glory. It is a powerful lesson: in computational science, how you measure can be just as important as what you measure.

So far, we have discussed systems where the future depends only on the present. The rate of change $\dot{\mathbf{x}}(t)$ is a function of the state $\mathbf{x}(t)$ at the very same instant. But many systems in nature, engineering, and economics have memory. Their evolution depends not only on their current state but also on their state at some time in the past. These are described by ​​Delay Differential Equations (DDEs)​​.

Think of a population of rabbits. The birth rate today doesn't depend on the number of rabbits today, but on the number of rabbits a gestation period ago. In control theory, the time it takes for a signal to travel from a sensor to an actuator introduces a delay. To compute the derivative at time $t$, the solver needs to know the solution at a past time, $t-\tau$, where $\tau$ is the delay.

This poses a huge challenge for a standard ODE solver. As it steps forward, say from $t_k$ to $t_{k+1}$, it needs the value of the solution at a past point $t-\tau$. But what if that point falls between the stored steps, for instance, between $t_{j}$ and $t_{j+1}$ from some earlier part of the integration? The solver needs to look back in its history and find the solution at a time it never explicitly calculated.

Dense output is the perfect "memory bank" for this task. By storing the polynomials that describe the solution over past intervals, the solver can accurately and efficiently "reconstruct" the solution at any past time it needs. It provides a continuous history of the trajectory, making the solution of DDEs feasible. The problem becomes even more fascinating with state-dependent delays, where the delay itself changes with the system's state, as in a hypothetical system like $y'(t) = -k \cdot y(t - c|y(t)|)$. Such problems are almost unimaginable to solve numerically without the continuous history provided by dense output.

We live in an increasingly digital world, where continuous physical processes are often monitored and controlled by discrete-time computers. Dense output plays a crucial role as a bridge between these two realms.

Consider a digital controller that samples the output of a continuous chemical plant every $T$ seconds. To design a high-performance controller, it's often not enough to know the state at the sampling instants. You might need to know the integral or the average value of a signal over the entire sampling period. A naive approach would be to assume the signal was constant or varied linearly between samples, but this introduces errors. The true behavior is a curve. As explored in advanced control problems, designing an exact digital equivalent of a continuous process requires knowing precisely how the continuous-time system behaves between the samples. For a general system, dense output provides the curve, which can then be integrated or processed to yield highly accurate digital equivalents.

This idea of "filling in the blanks" between samples creates a beautiful connection to a seemingly unrelated field: digital signal processing (DSP). In DSP, one often needs to increase the sampling rate of a signal—a process called interpolation. This is commonly done by inserting zeros between the original samples and then applying a digital low-pass filter to "fill in" the values of the newly created points. The process of designing this filter is a deep and fascinating subject.

What is the "best" way to fill in the blanks? In DSP, the answer is given by the Shannon-Nyquist sampling theorem. If the original signal was properly bandlimited, the ideal interpolation filter is a sinc function, $h[n] = \frac{\sin(\pi n/L)}{\pi n/L}$. This filter corresponds to a perfect "brick-wall" low-pass filter in the frequency domain, which flawlessly removes the spectral copies created by the upsampling process.

At first glance, this seems worlds away from our ODEs. But look closer. We have two fields, both trying to reconstruct a continuous reality from discrete points.

• The ​​ODE solver​​ uses dense output, often a piecewise polynomial like a cubic spline, to interpolate. Its "knowledge" comes from the dynamics of the system—the local derivatives—which dictate the shape of the curve.
• The ​​DSP system​​ uses a sinc-like filter. Its "knowledge" comes from the frequency content of the signal—the assumption that no high frequencies were present in the original signal.

Both are sophisticated forms of interpolation, each tailored to the fundamental principles of its domain. This parallel illustrates a profound unity in scientific computing: the challenge of representing the continuous with the discrete is universal, and the solutions, though born of different theories, share a common soul.

So, the next time you see the smooth curve of a simulated trajectory, remember the silent, sophisticated machinery of dense output working between the points. It is the artist that paints the full picture, the storyteller that fills in the narrative between the major plot points, revealing that in science, as in life, the real richness is often found not at the destinations, but in the journey between them.