Why Discontinuities Need Special Treatment
A delay differential equation (DDE) can carry information from an earlier time into its current rate of change. Consequently, an abrupt change in an input or a mismatch in the initial history can affect the solution again after a delay. A numerical solver that takes a step across one of these times may approximate a smooth curve where the true solution has a change in derivative. The resulting error can persist in later steps, including through the delayed values used by the model.
The aim here is to identify those times and make them boundaries of the numerical integration. The discussion concerns retarded DDEs with fixed positive delays, whose right-hand sides depend on the current state and delayed states, but not on delayed derivatives. Assume the right-hand side is smooth between known input events, and that the history is piecewise smooth. These conditions allow the solution to be treated as smooth between the event times and their delayed propagations.
For example, consider a system with fixed delays \(\tau_j>0\):
Here \(x(t)\) may be a vector of states, \(u(t)\) is an input, and a history function specifies \(x(t)\) for times before the start of integration. An abrupt input change can make \(x'\) discontinuous immediately. A discontinuity in \(x'\) can then affect \(x''\) when the corresponding delayed state is encountered. For smooth \(f\), a discontinuity in a derivative of order \(r\) can therefore produce a discontinuity in the derivative of order \(r+1\) one delay later. This is a propagation rule for candidate times; a particular effect can vanish if the relevant sensitivity is zero or different contributions cancel.
Propagation and Integration Boundaries
Suppose there is a single delay \(\tau\), and \(x'\) has a jump at time \(s\). At time \(s+\tau\), the delayed argument \(x(t-\tau)\) reaches \(x(s)\). The state itself is continuous at \(s\), so this does not generally make \(x'\) jump again. But the derivative of the delayed state, \(x'(t-\tau)\), has a jump there. Differentiating the DDE within smooth intervals shows how that jump enters \(x''\).
In this expression, \(f_z\) is the derivative of \(f\) with respect to its delayed-state argument, and \([g]_s=g(s+)-g(s-)\) denotes the right-hand limit minus the left-hand limit. The jump relation assumes the other terms are continuous at \(s+\tau\), with no simultaneous event contributing another jump. With multiple delays, each delayed-state derivative contributes a corresponding term. The relation explains both the one-order increase in derivative order and why some candidate discontinuities may be absent in a particular model.
The integration plan follows from these candidate times. Start with known events: input jumps, the initial time, and any known discontinuities in the history relevant to the integration interval. For each candidate time \(d\), add \(d+\tau_j\) for every delay, retaining times within the simulation interval. Repeat until no more times are needed. With multiple fixed delays, this creates sums of delays from the original event times. If different sums coincide, the time is included once, though the separate effects at that time still contribute to the net behavior.
Record input events and times where the history or initial conditions introduce a loss of smoothness.
Add each fixed delay to each known candidate time, repeatedly, up to the simulation end time.
End a solver segment at each candidate time and restart on the other side using the correct history values.
Inspect behavior near boundaries and repeat the simulation with tighter tolerances or smaller maximum steps.
A solver should not use a single smooth interpolant across a known discontinuity. Store or construct the history separately on each smooth segment, and make sure evaluations at a boundary use the appropriate one-sided behavior where needed. For a continuous state with a derivative jump, the state value at the boundary agrees from both sides, but a polynomial interpolant spanning the boundary can still distort nearby derivatives and delayed evaluations.
Worked Example: A History Mismatch and Its Propagation
Consider the hypothetical scalar DDE
Time is measured in days and \(x\) is dimensionless, so both coefficients have units of day\(^{-1}\). The history is constant, with left derivative zero at time zero. The DDE gives the right derivative
Thus \(x\) is continuous at zero, but \(x'\) has a jump of \(-0.1\) day\(^{-1}\). On the first interval, \(0<t<2\), the delayed value is the constant history \(1\). Solving the resulting ordinary differential equation gives
At \(t=2\), the delayed state \(x(t-2)\) reaches the initial point. The state and \(x'\) remain continuous there because the delayed state itself is continuous. Differentiating within the intervals gives \(x''=-0.4x'+0.3x'(t-2)\). Therefore
The candidate times are \(0,2,4,6,\ldots\) days. The solver should end a segment at each of these times within the simulation horizon. The derivative order affected can increase at each propagation, even though the state itself remains continuous.
Worked Example: Propagation from an Input Step
Suppose a hypothetical amount \(x\), measured in arbitrary amount units, follows
with zero history. Time is in days; \(u\) has units of amount per day, and the coefficients have units of day\(^{-1}\). Before day 1 the state is zero. The input step makes \(x'\) jump by \(2\) amount units per day at day 1. For \(1<t<2\), the delayed state is still zero, so
At day 2, \(x(t-1)\) reaches \(x(1)\), which is continuous. Hence the input step has not caused a jump in \(x'\) at day 2. But differentiating the equation gives \(x''=-0.2x'+0.1x'(t-1)+u'(t)\) between events. The delayed derivative carries the jump from day 1:
The input event is at day 1, and its successive propagated candidates occur at days 2, 3, and later. A suitable integration plan splits at these times, as well as at any other input events. Integrating across day 1 without representing the input change as an event would also give the wrong forcing over part of a solver step.
Worked Example: Candidate Times With Two Delays
Suppose a model has fixed delays of \(0.75\) day and \(1.25\) days, and an initial derivative mismatch at day 0. To plan a simulation through day 3, begin with the initiating time 0 and repeatedly add either delay. The candidate times are the nonnegative sums \(0.75m+1.25n\) that do not exceed 3, where \(m\) and \(n\) are nonnegative integers.
For \(n=0\), the times are 0, 0.75, 1.50, 2.25, and 3.00. For \(n=1\), they are 1.25, 2.00, and 2.75. For \(n=2\), the only additional time not exceeding 3 is 2.50. After sorting and removing duplicates, the integration boundaries are
The time 2.00 days has two routes: \(0.75+1.25\) and \(1.25+0.75\). It is one solver boundary, but effects arriving by both routes can add or cancel. These are candidate times, not a guarantee that a visible jump occurs at every one: the model’s delayed-state sensitivities and the incoming jumps determine whether each candidate produces a nonzero discontinuity.
Choosing and Checking a Numerical Strategy
An event-aware DDE solver or a method-of-steps implementation can both handle these problems, provided they represent discontinuities correctly. The method of steps solves ordinary differential equations on successive intervals while supplying already-computed delayed values. Its intervals must be split further when a propagated discontinuity lies inside the next interval. A general adaptive ODE solver with dense output may also be used, but it needs the same event schedule and piecewise history handling. Adaptively reducing the step size is helpful for accuracy; it does not, by itself, ensure that a step ends exactly at a known nonsmooth time.
| Approach | Strength | Risk or limitation |
|---|---|---|
| Event-aware segmentation | Respects known changes in smoothness and keeps interpolation within smooth intervals. | Requires a correct list of event and propagation times. |
| Adaptive steps without event boundaries | Can reduce local error in smooth regions and is convenient for general systems. | May step across a kink; the error estimate can miss the effect on delayed interpolation. |
| Uniformly very small steps | Simple to implement and useful as a rough check. | Can be inefficient and still does not align steps exactly with events. |
Check the schedule against every input discontinuity and every relevant delay. Around a candidate time, compare left- and right-side behavior at the expected derivative order rather than expecting the state itself to jump. Then compare solutions using tighter tolerances and smaller maximum step sizes. A stable state trajectory together with a changing derivative near a boundary may indicate unresolved nonsmoothness rather than a biological effect. If results shift materially under refinement, or solver warnings appear near boundaries, revisit both the event schedule and the representation of delayed history.
The propagation rule relies on a fixed-delay retarded DDE and smooth dynamics between events. A model with state resets can have actual state jumps; its reset rules must be included explicitly. A model with state-dependent delays can have event times that must be located as part of the computation. Delayed derivatives, as in neutral DDEs, require different regularity considerations. In all cases, a fine grid alone is not a substitute for representing known discontinuities.
Common Pitfall
A common mistake is to include an input jump at time \(t_0\) but omit its delayed consequences, reasoning that the state is continuous and therefore the later solution must be smooth. Continuity of \(x\) does not imply continuity of \(x'\), and continuity of \(x'\) at a propagated time does not imply continuity of \(x''\). The correction is to track the derivative order of the original event and propagate candidate times by adding each delay, while checking which derivative can be affected. Conversely, do not assert that every candidate time has a nonzero jump: sensitivities may be zero or effects may cancel.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- A scalar DDE has a delay of 1.5 days and an input jump at day 2. Give the first three candidate times, including the input event, at which the jump or its delayed effects may matter.
- In a retarded DDE, a derivative jump occurs at time \(s\). Under smoothness assumptions and absent coincident events, which derivative can acquire a jump at \(s+\tau\)?
- A model has delays of 1 day and 2 days, with an initiating event at time 0. List the distinct candidate times through day 4.
- Why can a solver that takes a very small step still be inaccurate near a known discontinuity?
Solutions. (1) The event is at day 2; adding the 1.5-day delay repeatedly gives candidate times 3.5 and 5.0 days as the next two propagated times. (2) The second derivative \(x''\) can acquire a jump, because the delayed derivative enters the differentiated DDE. (3) The distinct sums through day 4 are 0, 1, 2, 3, and 4 days. (4) A small step may still cross the nonsmooth time, causing the solver or its dense interpolant to approximate behavior on both sides as though it were smooth. Explicitly ending a segment at the event avoids that particular error.