Streamline Integration in Vector Fields and Where the Error Accumulates

A vector field assigns a direction to every point, and the natural picture of it is a set of curves tangent to those directions everywhere. The curves are defined by an ordinary differential equation, $\mathrm{d}\mathbf{x}/\mathrm{d}s = \mathbf{v}(\mathbf{x})$, integrated from a starting point, and a reader looking at the result is implicitly making a statement about the field’s behaviour over the whole interval the curve covers. A line that drifts off a real streamline is not a slightly wrong line; it shows a trajectory that does not exist.

Two sources of error produce that drift, and they are quite different in character. The integrator has a truncation error that shrinks predictably with step size. The field lookup has an interpolation error that does not shrink with step size at all, because refining the integration samples the same reconstructed field more precisely. For most datasets plotted at ordinary resolution, the second dominates.

The field that is actually integrated is a reconstruction

A vector field stored on a grid exists as samples. The integrator queries it at arbitrary positions, so it evaluates the interpolation model, and the streamline follows the integral curves of that model, not of the sampled field and certainly not of the underlying continuous field. This distinction explains a family of behaviours that otherwise look like bugs.

Inside a cell, the reconstructed field has a specific form. On a simplicial grid, a linear interpolation on each triangle or tetrahedron makes the velocity components affine functions of position, so the trajectory within a cell is the orbit of a linear system, $\dot{\mathbf{x}} = \mathbf{A}\mathbf{x} + \mathbf{b}$, whose exact solution is a matrix exponential. For a two-dimensional affine field the orbits are conic sections: ellipses when the matrix has complex eigenvalues, hyperbolic arcs when the eigenvalues are real and of opposite sign, and straight lines in the degenerate cases. A streamline crossing such a cell is therefore an arc of a conic, not an arbitrary smooth curve. On a hexahedral grid the reconstruction is trilinear rather than affine, and the trajectory has a more complicated closed form, but the qualitative consequence is the same: the curve being approximated belongs to the interpolation model.

At a cell boundary, the reconstructed field is continuous in value but its derivative jumps. The integrator’s error estimator sees this as a change in the local solution behaviour, and adaptive controllers react by shrinking the step exactly where the field is least smooth. Worse, the reconstruction’s interpolation error makes the reconstructed field’s critical points (where $\mathbf{v} = 0$) sit at positions determined by the grid, not by the data. A saddle point in the real field is displaced, and near it streamlines are highly sensitive: a small displacement of the critical point changes which quadrants streamlines enter and leave, and therefore changes the topology of the picture.

Step size, order, and the saturation point

A fourth-order Runge-Kutta step has a local truncation error of order $h^5$ and a global error of order $h^4$. Halving the step divides the global error by roughly sixteen for a smooth field. The error in the streamline’s position after integrating over a length $L$ therefore improves rapidly, and the obvious response to a visibly polygonal streamline is to reduce $h$.

That response stops paying once the interpolation error dominates. The interpolation error of a trilinear reconstruction is second order in the cell spacing and is a property of the data, not of the integration. Once the integration error falls below it, further refinement traces the same wrong curve more precisely. The practical diagnostic is to integrate the same streamline with two step sizes and compare the resulting positions: if halving $h$ does not halve the separation, the limit is the field model.

Adaptive control is standard for this reason. A step-doubling or embedded pair estimates the local error, accepts the step if it is below tolerance, and otherwise retries with a smaller step. The step size becomes small in regions of high curvature and large in straight regions, which is the right allocation of work. The controller also needs an absolute floor: as the field magnitude goes to zero near a critical point, the parameter step needed to advance a fixed geometric distance grows without bound, and without a term that scales the tolerance by the local speed, the integration slows to a crawl or the step grows to the point where the streamline exits the cell entirely.

A related trap is normalization. Integrating $\mathrm{d}\mathbf{x}/\mathrm{d}s = \mathbf{v}/\lVert\mathbf{v}\rVert$ makes the parameter an arc length and gives a fixed number of steps per unit distance, which is convenient for a fixed-length streamline. It also removes the magnitude information from the parameterization, so a streamline through a slow region covers the same geometric distance per step as one through a fast region, and the step count no longer reflects anything about the field. For visualization purposes that is usually desirable; for analysis, where the travel time matters, it destroys the quantity of interest.

import numpy as np

def rk4_streamline(field, x0, h, max_steps, bounds=None, h_min=1e-6):
    """Integrate dx/ds = field(x) with a fixed step, stopping on bounds or stagnation."""
    pts = [np.asarray(x0, dtype=float)]
    x = pts[0].copy()
    for _ in range(max_steps):
        k1 = field(x)
        speed = np.linalg.norm(k1)
        if speed == 0.0:
            break
        step = min(h, 0.5 / speed)         # cap geometric advance per step
        step = max(step, h_min)
        k2 = field(x + 0.5 * step * k1)
        k3 = field(x + 0.5 * step * k2)
        k4 = field(x + step * k3)
        x = x + (step / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
        if bounds is not None and not bounds.contains(x):
            break
        pts.append(x.copy())
    return np.array(pts)

The cap on geometric advance is the coarse version of a curvature-aware controller. It prevents a step from leaping over a cell boundary in a single jump, which is the specific failure that makes a streamline appear to tunnel through a region without following it.

Seeding is a sampling problem, not a decoration

Given a field and an integrator, the picture is still undetermined, because the set of streamlines drawn is chosen by the seed positions. Uniform seeding over a plane produces dense lines where the field converges and sparse lines where it diverges, since lines seeded in a converging region crowd together downstream. Distance-based seeding places new seeds only where the existing lines are farther apart than a threshold, measured in the image plane or in a plane perpendicular to the flow. The threshold translates directly into a density and therefore into the visual load.

The trade is between coverage and clutter. A threshold that is too loose leaves regions unrepresented and a reader infers a featureless field where the sampling simply missed; one that is too tight obscures the underlying structure with ink. Since the density of lines in the image scales as the inverse square of the spacing threshold in three dimensions, doubling the apparent detail quadruples the line count and the integration cost.

In three dimensions, the corresponding object is a stream surface, swept by integrating a curve of seeds and connecting the resulting streamlines into a ribbon. The surface has the advantage of occluding less than a solid and revealing twisting, but it can self-intersect, which has no physical meaning for a steady flow and is purely an artefact of the projection. Illumination along the surface then becomes the cue that carries the twist.

Integration direction matters as well. A steady field defines a trajectory forward and backward from a seed, and a picture that shows only forward lines truncates the structure near the seed plane, hiding the upstream part of the same feature. Both directions are usually integrated, with the backward result reversed and concatenated.

Termination, and the cases where the equation is not defined

Every streamline integrator needs a stop condition, and each of the standard ones corresponds to a physical or numerical boundary.

Leaving the domain is the simple case: once the position exits the data bounds, the field has no definition there, and continuing requires extrapolation that the data does not support. The clamp-versus-missing choice from Sample Positions, Trilinear Interpolation, and the Half-Cell Offset applies directly, and clamping here is dangerous because it creates a boundary layer of artificial flow along the domain edge that draws streamlines parallel to the wall.

Stagnation is where the speed falls below a threshold. In a field with a true critical point, the trajectory approaches it asymptotically and never arrives, so without a threshold the integration runs until the step cap. The threshold converts an infinite-time approach into a finite stop, and its value determines how far the line is drawn toward the point. Too large and the line stops short, leaving a visible gap around every critical point; too small and the line fades into the numerical noise near the point, where the interpolation error is comparable to the field’s magnitude.

Recirculation is the third case: a closed orbit, a vortex core, or a periodic domain. A fixed step count ends the line after an arbitrary number of turns, and a line count based on arc length ends it at an arbitrary phase. Detecting closure requires comparing the current position and direction against the seed history with a tolerance, and in a field integrated with interpolation error the orbit does not close exactly, so the criterion must be a band rather than a point. Failing to detect it produces a line that spirals slowly across the plane instead of closing, which for a periodic domain looks like a drift artefact.

Where the field is not steady, the equation gains a time dependence and the picture changes kind. A path line follows one particle over time and requires the field at each instant; a streak line connects particles released earlier and is what a dye injection shows. Neither is a streamline of the instantaneous field, and the three coincide only for a steady flow. A visualization that labels a path line as a streamline is making a claim the data does not support, and in an unsteady field the difference is qualitative, not cosmetic.

Where the model runs out

Interpolation-based streamlines inherit every limitation of the reconstruction. A field sampled at a resolution coarser than the structures it contains cannot be fixed by better integration; the structures are not in the data. Adaptive mesh refinement helps only where the mesh was refined, and a streamline crossing from a refined region into a coarse one changes its error character mid-path.

Vector field topology, the skeleton of critical points and their separatrices, is the compact summary that does not depend on seeding, but it depends even more sharply on the interpolation model, since critical point positions are exactly the places where the reconstruction error is largest in relative terms. Computing topology from a reconstruction is therefore a statement about the reconstruction.

The scalar counterpart of the same sampling problem is volume rendering, where a ray also crosses cells and the reconstruction determines what the ray accumulates; the transfer function that turns those samples into color and opacity is the subject of Direct Volume Rendering and the Opacity Correction Nobody Applies. When the field is produced by a running simulation rather than stored, the integration has to happen while the data is still resident, which is the constraint described in In-Situ Analysis and the Bandwidth That Is Not There.