Handling Discontinuities and Piecewise Correlations in an Engineering Solver

Many engineering models do not behave smoothly from start to finish. A valve opens, a controller saturates, a material property changes phase, or an empirical correlation switches from its laminar form to its turbulent form. From the point of view of the numerical solver, these are discontinuities: points where the governing equations, or the solution itself, change abruptly. Because most integration algorithms reconstruct the solution by fitting polynomials through earlier values, an unhandled break in smoothness can degrade accuracy, shrink the step size, or stall convergence altogether. Solvers deal with the problem using two complementary strategies: treating the switch as an event that is located precisely and used to restart the integration, or replacing the sharp transition with a smoothed function. This article explains how both approaches work, how they apply to piecewise-defined correlations, and how to troubleshoot the convergence problems that discontinuities cause.

What Counts as a Discontinuity in a Solver?

In an ordinary differential equation (ODE) solver, a discontinuity is any point where the smoothness assumptions behind the numerical method break down. Engineering models introduce them constantly. Cellier’s lecture notes on continuous system modeling list switching, limiters, dry friction, and impulse-like phenomena as typical sources, and emphasize that these features strongly influence the numerical behavior of the underlying differential equation solver.

Not all discontinuities are the same, and the distinction matters because each type calls for different treatment. The documentation of the AMICI simulation package, which handles events in dynamical systems, distinguishes three types:

• Solution jump discontinuities. The solution itself jumps. Examples include SBML events (the Systems Biology Markup Language’s mechanism for resetting model quantities at trigger conditions) or delta functions in the right-hand side of the equations.
• Right-hand-side (RHS) jump discontinuities. The equations change their branch, but the solution remains continuous — hence the discontinuity is removable in the solution. These are created by piecewise, Heaviside, and other logical operators in the right-hand side.
• RHS removable discontinuities. Neither the solution nor its first derivative jumps, but higher-order temporal derivatives do. These arise from functions such as max or min used to clamp or select values in the equations.

Solution jumps require an explicit reset of the state when the event fires. The other two types are handled by changing which branch of the equations is active — but as discussed below, even a fully continuous solution with kinked derivatives can trouble an integrator.

Why Polynomial-Based Integrators Struggle

Most differential equation solvers currently on the market operate on polynomial extrapolation: at each step they fit a polynomial of some order through known values of the state variable and its derivative at the current and past times, then evaluate that polynomial one step ahead. The step-size and error controls are calibrated for smooth functions, where higher-order extrapolation is accurate and the error estimate is trustworthy.

A discontinuity inside a step invalidates that premise. If the model quietly switches branch while the integrator is straddling the switch, the polynomial is fitted across a break, the error estimate can be wrong, and the solver may either waste effort on rejected steps or advance past the switch without locating it. When parameter sensitivities are computed alongside the solution — as in parameter estimation — logical operations and events need dedicated formulas to keep those sensitivities correct, a point documented in the AMICI implementation notes and traced there to Fröhlich and colleagues’ work on parameter estimation with discrete events and logical operations.

Strategy 1: Detect the Event, Then Restart the Integration

The first strategy is to make the discontinuity explicit. A state event fires when a user-defined switching function — some quantity computed from the state and time — crosses a threshold, conventionally zero. A tank level reaching a limit, a switch closing at a computed current, or a friction law activating when relative velocity passes through zero are all state events.

The standard workflow is:

1. Integrate normally while monitoring each switching function for a sign change within the step.
2. When a sign change appears, locate the crossing time accurately by root-finding within the step.
3. Apply whatever the event demands: change the active branch of the equations, or jump the state itself.
4. Reinitialize and restart the integrator with fresh history from the event time.

Dedicated codes automate this. DISODE45 (ACM TOMS Algorithm 968, published 2016 by Calvo, Montijano, and Rández) is an adaptive Runge–Kutta code built on the DOPRI5(4) pair for piecewise-smooth initial value problems: it automatically detects and accurately locates switching points, restarting the integration after each discontinuity. The user supplies the functions defining the hypersurfaces on which switching points lie. Importantly, DISODE45 also handles Filippov systems — systems whose right-hand side differs on either side of a switching surface, and whose solutions may slide along that surface rather than crossing it. Dry friction is the classic engineering case: near stiction, a naive implementation can chatter between branches, while a sliding-mode treatment integrates a combined regime instead.

The practical lessons are to expose switching conditions as explicit functions the solver can root-find, to restart cleanly after every located event, and to apply any state reset before the restart.

Strategy 2: Smooth the Sharp Transition

The second strategy is to remove the mathematical sharpness where the physics allows it. COMSOL’s reference on specifying discontinuous functions notes that a step in a coefficient or material property causes real numerical trouble: in time-dependent studies the time-stepping algorithm can run into problems, and in stationary studies mesh-resolution issues arise, such as overshooting or undershooting of the solution near sharp fronts.

The conventional way to write a sharp transition is with a logical expression that evaluates to 1 when true and 0 otherwise — COMSOL’s documentation uses the example sin(2πt)·(t < 10), a sine wave that exists for ten seconds and is zero afterwards. COMSOL recommends replacing such discontinuities with a smoothed step function, for two reasons: numerical reliability and convergence improve, and the result is often more physical, because a step function that appears in a model is, in reality, a smoothed transition shaped by inertia.

Practically, this means that a limiter, saturation, or property switch should be represented by a function that ramp steeply but continuously across the transition, using the predefined Step function (which includes smoothing by default in COMSOL) or an equivalent construct in your tool. The caveat is that smoothing must remain narrow relative to the time scales of interest: too generous a smoothing width changes the model’s answers, while too narrow a width brings the numerical problems back.

Making Piecewise Correlations Solver-Friendly

A piecewise correlation is an empirical relationship supplied as different formulas over disjoint ranges of an independent variable or state. They are ubiquitous in flow and heat-transfer calculations. Pipe-flow friction factors, for example, take laminar forms at low Reynolds numbers and turbulent forms (such as the Blasius or Colebrook-type relations) at high Reynolds numbers, with the flow regime conventionally considered to change near a Reynolds number of about 2,300. Internal-flow heat-transfer correlations similarly change shape between laminar and turbulent regimes, and thermodynamic property formulations switch branch at saturation boundaries.

Piecewise correlations cause solver trouble when their branch boundaries are handled implicitly. Two failure modes are common: the branches do not agree at the switch (so the correlation itself jumps), or they agree in value but not in slope (leaving a kink that degrades high-order integration). The following practices keep them manageable:

• State the switching criterion explicitly and compute it as a smooth switching function the solver can root-find, rather than burying the decision inside conditional branches.
• Where the data permit, join branches so the correlation is continuous — ideally in both value and first derivative — at the regime boundary.
• If literature branches genuinely disagree at the switch, either accept the jump and handle it as an event, or blend the branches over a narrow band and document the band.
• Avoid implementing regime selection with nested max/min or Heaviside operators when sensitivities matter; AMICI’s documentation notes that these create removable discontinuities in higher derivatives and require special event-aware treatment.
• Document which branch applies exactly at the switch value (interval conventions), so different users and software versions get identical results.

A Practical Workflow for Solver Users

1. Inventory the conditional logic. List every switch, limiter, clamp, piecewise correlation, and event in the model.
2. Classify each item. Decide whether it is a solution jump, an RHS branch change, or only a kink in higher derivatives, following the three-type taxonomy above.
3. Expose switching functions. For state-dependent switches, define explicit functions of the state whose zeros the solver can locate.
4. Restart cleanly. After each located event, apply any state jump and reinitialize the integrator’s history before continuing.
5. Smooth deliberately. Apply smoothing only to transitions where no jump is modeled, keep the width narrow relative to the dynamics, and justify it physically.
6. Watch the diagnostics. Step-size histories and convergence logs concentrated near switch times usually reveal unhandled discontinuities.
7. Validate. Compare against smoothed and event-handled variants of the same model, and check sensitivities against finite-difference values if the model is used for parameter estimation.

Symptoms, Causes, and Fixes

Trade-Offs to Keep in Mind

Smoothing and event handling are complements, not rivals, but each has limits. Smoothing preserves a single continuous solve, yet it changes the model: too wide a smoothing band rounds away genuinely sharp behavior, and too narrow a band recreates the convergence problems. Event handling is faithful to the discontinuity but adds implementation effort and cost, and systems with dense sequences of events slow simulations; sliding modes additionally require specialized algorithms, as the DISODE45 design shows. Knowing which type of discontinuity you face is therefore the decisive first step, because it dictates which treatment is appropriate.

Key Takeaways

• Discontinuities from switching, limiters, dry friction, and impulses are normal in engineering models, but they strongly affect solver numerics.
• Classify each discontinuity: a jump in the solution needs an event-driven state reset; branch changes need located events or smoothing; clamps like max/min mostly hurt via higher derivatives and sensitivities.
• For piecewise correlations, make switch criteria explicit, join branches as smoothly as the data allow, and document conventions at the boundary.
• Use smoothed steps where inertia justifies them and keep their width narrow; use event detection with clean restarts where the switch is real.