Solving Delay Differential Equations

Delay differential equations (DDEs) are ordinary differential equations that relate the solution at the current time to the solution at past times. This delay can be constant, time-dependent, state-dependent, or derivative-dependent. In order for the integration to begin, you generally must provide a solution history so that the solution is accessible to the solver for times before the initial integration point.

Constant Delay DDEs

A system of differential equations with constant delays has the form:

$$y^{\prime }\left(t\right)=f\left(t,y\left(t\right),y\left(t-{\tau }_1\right),\dots ,y\left(t-{\tau }_k\right)\right).$$

Here, t is the independent variable, y is a column vector of dependent variables, and y ′ represents the first derivative of y with respect to t. The delays, τ₁,…,τ_k, are positive constants.

The dde23 function solves DDEs with constant delays with history y(t) = S(t) for t <t₀.

The solutions of DDEs are generally continuous, but they have discontinuities in their derivatives. The dde23 function tracks discontinuities in low-order derivatives. It integrates the differential equations with the same explicit Runge-Kutta (2,3) pair and interpolant used by ode23. The Runge-Kutta formulas are implicit for step sizes bigger than the delays. When y(t) is smooth enough to justify steps this big, the implicit formulas are evaluated by a predictor-corrector iteration.

Time-Dependent and State-Dependent DDEs

Constant time delays are a special case of the more general DDE form:

$$y^{\prime }\left(t\right)=f\left(t,y\left(t\right),y\left(d{y}_1\right),\mathrm{...},y\left(d{y}_p\right)\right).$$

Time-dependent and state-dependent DDEs involve delays dy₁,..., dy_k that can depend on both time t and state y. The delays dy_j(t, y) must satisfy dy_j(t, y) ≤ t on the interval [t₀, t_f] with t₀ < t_f.

The ddesd function finds the solution, y(t), for time-dependent and state-dependent DDEs with history y(t) = S(t) for t < t₀. The ddesd function integrates with the classic four-stage, fourth-order explicit Runge-Kutta method, and it controls the size of the residual of a natural interpolant. It uses iteration to take steps that are longer than the delays.

DDEs of Neutral Type

Delay differential equations of neutral type involve delays in y ′ as well as y:

$$y^{\prime }\left(t\right)=f\left(t,y\left(t\right),y\left(d{y}_1\right),\mathrm{...},y\left(d{y}_p\right),y^{\prime }\left(dy{p}_1\right),\mathrm{...},y^{\prime }\left(dy{p}_q\right)\right).$$

The delays in the solution must satisfy dy_i(t,y) ≤ t. The delays in the first derivative must satisfy dyp_j(t,y) < t so that y ′ does not appear on both sides of the equation.

The ddensd function solves DDEs of neutral type by approximating them with DDEs of the form given for time-dependent and state-dependent delays:

$$y^{\prime }\left(t\right)=f\left(t,y\left(t\right),y\left(d{y}_1\right),\mathrm{...},y\left(d{y}_p\right)\right).$$

For more information, see Shampine [1].

Evaluating the Solution at Specific Points

Use the deval function and the output from any of the DDE solvers to evaluate the solution at specific points in the interval of integration. For example, y = deval(sol, 0.5*(sol.x(1) + sol.x(end))) evaluates the solution at the midpoint of the interval of integration.

History and Initial Values

When you solve a DDE, you approximate the solution on an interval [t₀,t_f] with t₀ < t_f. The DDEs show how y ′(t) depends on values of the solution (and possibly its derivative) at times prior to t. For example, with constant delays y ′(t₀) depends on y(t₀ – τ₁),…,y(t₀ – τ_k) for positive constants τ_j. Because of this, a solution on [t₀, t_k] depends on values it has at t ≤ t₀. You must define these values with a history function, y(t) = S(t) for t <t₀.

Discontinuities in DDEs

If your problem has discontinuities, it is best to communicate them to the solver using an options structure. To do this, use the ddeset function to create an options structure containing the discontinuities in your problem.

There are three properties in the options structure that you can use to specify discontinuities; InitialY, Jumps, and Events. The property you choose depends on the location and nature of the discontinuities.

Nature of Discontinuity	Property	Comments
At the initial value t = t0	InitialY	Generally the initial value y(t0) is the value S(t0) returned by the history function, meaning the solution is continuous at the initial point. If this is not the case, supply a different initial value using the InitialY property.
In the history, i.e., the solution at t <t0, or in the equation coefficients for t >t0	Jumps	Provide the known locations t of the discontinuities in a vector as the value of the Jumps property. Applies only to dde23.
State-dependent	Events	dde23, ddesd, and ddensd use the events function you supply to locate these discontinuities. When the solver finds such a discontinuity, restart the integration to continue. Specify the solution structure for the current integration as the history for the new integration. The solver extends each element of the solution structure after each restart so that the final structure provides the solution for the whole interval of integration. If the new problem involves a change in the solution, use the InitialY property to specify the initial value for the new integration.

Propagation of Discontinuities

Generally, the first derivative of the solution has a jump at the initial point. This is because the first derivative of the history function, S(t), generally does not satisfy the DDE at this point. A discontinuity in any derivative of y(t) propagates into the future at spacings of τ₁,…, τ_k when the delays are constant. If the delays are not constant, the propagation of discontinuities is more complicated. For neutral DDEs of the forms in Constant Delay DDEs and Time-Dependent and State-Dependent DDEs, the discontinuity appears in the next higher order derivative each time it is propagated. In this sense, the solution gets smoother as the integration proceeds. Solutions of neutral DDEs of the form given in DDEs of Neutral Type are qualitatively different. The discontinuity in the solution does not propagate to a derivative of higher order. In particular, the typical jump in y ′(t) at t₀ propagates as jumps in y ′(t) throughout [t₀, t_f].

DDE Examples and Files

Several available example files serve as excellent starting points for most common DDE problems. To easily explore and run examples, simply use the Differential Equations Examples app. To run this app, type

odeexamples

To open an individual example file for editing, type

edit exampleFileName.m

To run an example, type

exampleFileName

This table contains a list of the available DDE example files, as well as the solvers and the options they use.

Example File	Solver Used	Options Specified	Description	Example Link
ddex1	dde23	—	DDE with constant history	DDE with Constant Delays
ddex2	dde23	'Jumps'	DDE with a discontinuity	Cardiovascular Model DDE with Discontinuities
ddex3	ddesd	—	DDE with state-dependent delays	DDE with State-Dependent Delays
ddex4	ddensd	—	Neutral DDE with two delays	DDE of Neutral Type
ddex5	ddensd	—	Neutral DDE with initial value	Initial Value DDE of Neutral Type

References

[1] Shampine, L.F. “Dissipative Approximations to Neutral DDEs.” Applied Mathematics & Computation, Vol. 203, 2008, pp. 641–648.