Delay-Differential Equations

• Dr. Skip Thompson, Radford University, Radford, Virginia

Delay differential equations differ from ordinary differential equations in that the derivative at any time depends on the solution (and in the case of neutral equations on the derivative) at prior times. The simplest constant delay equations have the form \[\tag{1} y'(t) = f(t, y(t), y(t-\tau_1), y(t-\tau_2),\ldots, y(t-\tau_k)) \]

where the time delays (lags) \( \tau_j \) are positive constants. More generally, state dependent delays may depend on the solution, that is \( \tau_i = \tau_i (t,y(t)) \ .\)

Introduction

Systems of delay differential equations now occupy a place of central importance in all areas of science and particularly in the biological sciences (e.g., population dynamics and epidemiology). Baker, Paul, & Willé (1995) contains references for several application areas.

Interest in such systems often arises when traditional pointwise modeling assumptions are replaced by more realistic distributed assumptions, for example, when the birth rate of predators is affected by prior levels of predators or prey rather than by only the current levels in a predator-prey model. The manner in which the properties of systems of delay differential equations differ from those of systems of ordinary differential equations has been and remains an active area of research; see Martin & Ruan (2001) and Raghothama & Narayanan (2002) for typical examples of such studies. See also Shampine, Gladwell, and Thompson (2003) for a description of several common models.

Initial History Function

Additional information is required to specify a system of delay differential equations. Because the derivative in (1) depends on the solution at the previous time \( t - \tau_j \ ,\) it is necessary to provide an initial history function to specify the value of the solution before time \( t = 0 \ .\) In many common models the history is a constant vector; but non constant history functions are encountered routinely. For most problems there is a jump derivative discontinuity at the initial time. In addition, any solution or derivative discontinuity in the history function at points prior to the initial time need to be handled appropriately since such discontinuities are propagated to future times.

Derivative Discontinuities

In most models, the delay differential equation and the initial history are incompatible: for some derivative order, usually the first, the left and right derivatives are not equal. For example, the simple model \( y'(t) = y(t-1) \) with constant history \( y(t) = 1 \) has the property that \( y'(0^{+}) = 1 \ne y'(0^{-}) = 0 \ .\)

One of the most fascinating properties of delay differential equations is the manner in which such derivative discontinuities are propagated in time. For the equation and history just described, for example, the initial first discontinuity is propagated as a second degree discontinuity at time \( t = 1 \ ,\) as a third degree discontinuity at time \( t = 2 \ ,\) and, more generally, as a discontinuity in the \( {(n+1)}^{st}\) derivative at time \( t = n \ .\) This behavior is typical of that for a wide class of delay differential equations: generalized smoothing occurs as the initial derivative discontinuity is propagated successively to higher order derivatives. Smoothing need not occur for neutral equations or for non-neutral equations with vanishing delays.

Neves & Feldstein (1976) characterized the tree of derivative discontinuity times for state dependent delay differential equations as the zeroes with odd multiplicity of equations \[\tag{2} \tau_i (t,y(t)) - T = 0 \]

where \( T \) is the initial time or any later discontinuity time.

Continuous Extensions

Several of the solvers discussed in the next section use explicit Runge-Kutta methods to integrate systems of delay differential equations. An important question in this case is that of interpolation. Unlike ordinary differential equation solvers that are based on linear multistep methods possessing natural extensions, early Runge-Kutta solvers did not incorporate interpolation; rather they stepped exactly to the next output point instead of stepping beyond it and obtaining interpolated solutions. Interest in the issues of obtaining dense output without limiting the step size in this fashion and by the desire to incorporate root finding led to the development of Runge-Kutta methods endowed with suitable interpolants. Interpolation is handled in one of two ways in modern Runge-Kutta solvers, Hermite interpolation and continuously imbedded methods. For example, the solver dde23 which is based on a third order Runge-Kutta method uses Hermite interpolation of the old and new solution and derivative to obtain an accurate interpolant. By way of contrast, the solver dde_solver uses a sixth order Runge-Kutta method based on a continuously embedded \( C^1 \) interpolant derived from the same derivative approximations used by the basic method. In addition to providing accurate and efficient solutions, either type of interpolant can be used in conjunction with a root finder to locate derivative discontinuity times.

Available Delay Differential Equation Software

A number of issues must be taken into account by software for delay differential equations. Baker, Paul, & Willé (1995), Shampine & Thompson (2001), and Thompson & Shampine (2006) discuss the various issues. The well known dmrode solver (Neves (1975)) was the first effective software for delay differential equations. Many of the central ideas on which this solver was based were used in later f77 solvers dklag5 (Neves & Thompson (1992)) and dklag6 (Corwin, Sarafyan, and Thompson (1997)), and the Fortran 90/95 dde_solver (Thompson & Shampine (2006)). Although the state of the art for numerical software for delay differential equations is not as advanced as that for ordinary differential equation software, several high quality solvers have recently been developed. The effectiveness of the software is determined in large part by the manner in which propagated derivative discontinuities are handled. Some delay differential equation solvers such as those in Paul (1995), and Thompson & Shampine (2006) explicitly track and locate the zeroes of (2) and include them as integration mesh points. Different approaches are used in other software. For example, the ddverk solver (Enright & Hayashi (1997)) uses derivative defect error control to implicitly locate discontinuity times. It then uses special interpolants to step cross the discontinuities. The ddesd solver (Shampine (2005)) uses residual error control to avoid the use of embedded local error estimates near discontinuity times.

Effective delay differential equation software must deal with other difficulties peculiar to systems of delay differential equations. Early software, for example, limited the step sizes used to be no larger than the smallest delay. But small delays are encountered in many problems; and this artificial restriction on the step size can have a drastic effect on the efficiency of a solver. Most of the solvers mentioned above are based on pairs of explicit continuously embedded Runge-Kutta methods (Shampine (1994)). When the step size exceeds a delay, the underlying interpolation polynomials are iterated in a manner somewhat akin to a predictor-corrector iteration for linear multistep methods. Refer to Baker & Paul (1996), Baker, Paul, & Willé (1995), Enright & Hayashi (1998), and Shampine & Thompson (2001) for details of various aspects of this issue.

The solvers dde23, ddesd, and dde_solver contain a very useful provision for finding zeroes of event functions (Shampine (1994)) that depend on the solution. In addition to solving a system of delay differential equations, they simultaneously locate zeroes of state dependent functions \( g(t,y(t)) = 0 \ .\) Such special events may signal problem changes requiring integration restarts. The use of event functions is illustrated in the next section.

Although much recent delay differential equation software utilizes explicit continuously embedded Runge-Kutta methods, software based on other methods has been developed. For example, Jackiewicz & Lo (2006) and Willé & Baker (1992) utilize generalized Adams linear multistep methods; and the radar5 solver (http://www.unige.ch/~hairer/software.html) is based on collocation methods. Another well known and widely used program with the ability to solve delay differential equations is the xppaut program (Ermentrout (2002)). The use of software based on a class of general linear methods (diagonally implicit multistage integration methods) is discussed in Hoppensteadt & Jackiewicz (2006) in conjunction with the problem considered in the next section. Bellen & Zennaro (2003) discuss the commonly used methods for delay differential equations in considerable detail.

An Example

Hoppensteadt & Jackiewicz (2006) investigated a model which generalizes previously studied models for infectious diseases. Solving this model requires the determination of a threshold time at which the accumulated dosage of infection reaches a prescribed level. Once this time is determined, the relevant equations may be integrated to obtain the desired solution. The minimum threshold time \( t_0 \) is the unique value for which \[ \int_{0}^{t_0} \rho(t) I_0(t) dt = m . \] The defining delay differential equations are \[ \tau'(t) = \frac{\rho(t) I(t)}{\rho(\tau(t)) I(\tau(t))}, \quad \tau(t_0) = 0 \] \[ S'(t) = -r(t) I(t) S(t), \quad S(0) = 0 \] Here the function \( I(t) \) is \[ I(t) = \begin{cases} I_0(t), & -\sigma \le t \le t_0, \\ I_0(t) + S_0 - S(\tau(t)), & t_0 \le t \le t_0 + \sigma \\ S(\tau(t-\sigma)) - S(\tau(t)), & t_0 + \sigma \le t \end{cases} \]

For Example 1 of the reference, the relevant variables and functions are given by \( m = 0.1, \sigma = 1, S_0 = 10, \rho(t) = 1, r(t) = r_0, S_0 = 10, \) and \[ I_0 = \begin{cases} 0.4(1+t), & t \le 0 \\ 0.4(1-t), & 0 \le t \le 1 \\ 0, & 1 < t \end{cases} \]

Following is a brief description of the manner in which this problem can be solved. As the description suggests, considerable dexterity may be required to solve a realistic system of delay differential equations. The solution of this problem involves three solution phases. Three delay differential equations are solved in each phase, one for \( \tau'(t) \ ,\) one for \( S'(t) \ ,\) and one for the accumulated dosage. The accumulated dosage \( y_3(t) \) is obtained by solving the equation \[ y_3'(t) = \rho(t) I_0(t) . \] Three additional delay functions, \( t-1 \ ,\) \( \tau(t) \ ,\) and \( \tau(t) - \sigma \ ,\) can be used to facilitate interpolations that must be performed during the different phases of the solution. During the first phase of the integration, \( \tau'(t) = 0 \ .\) A solution dependent event function \( g(t) = y_3(t) - m \) may be used to locate the threshold time \( t_0 \ ,\) i.e., the time \( t_0 \) for which \( y_3(t_0) = m \ .\) Once \( t_0 \) is located, the relevant equations for the second phase of the solution can be solved. The event function is changed to \( g(t) = t - (t_0 + \sigma) \ .\) This allows the program to locate \( t_0 + \sigma \) as an event. Once the second event is located, the full set of equations can be integrated to the desired final time of \( t = 8 \ .\) Two examples are considered in the paper. For the first example, \( \rho(t) = 1 \) and \( r(t) = r_0 \ .\) For the second example, \( \rho(t) = e^{-t} \) and \( r(t) = r_0 (1 + \sin(5t)) \ .\) Plots of \( S(t) \) and \( I(t) \) are given for several values of the parameter \( r_0 \ .\) To obtain these graphs, the problem is solved for each value of \( r_0 \ .\) Plots of the solutions thus obtained are displayed in Figures 1 and 2; they replicate those given in the paper.

Figure 1: Infected Individuals

Figure 2: Susceptible Individuals

References

• Baker, C.T.H. and Paul, C.A.H., A global convergence theorem for a class of parallel continuous explicit Runge-Kutta methods and vanishing lag delay differential equations, SIAM J. Numer. Anal., 33:1559-1576, 1996.

• Baker, C.T.H., Paul, C.A.H., and Willé, D.R., A bibliography on the numerical solution of delay differential equations, Numerical Analysis Report 269, Mathematics Department, University of Manchester, U.K., 1995.

• Baker, C.T.H., Paul, C.A.H., and Willé, D.R., Issues in the numerical solution of evolutionary delay differential equations, Adv. Comput. Math., 3:171-196, 1995.

• Bellen, A. and Zennaro, M., Numerical Methods for Delay Differential Equations, Oxford Science Publications, Clarendon Press, Oxford, 2003.

• Corwin, S.P., Sarafyan, D. and Thompson S., DKLAG6: a code based on continuously imbedded sixth order Runge-Kutta methods for the solution of state dependent functional differential equations, Appl. Numer. Math., 24:319-333, 1997.

• Enright, W.H. and Hayashi, H., A delay differential equation solver based on a continuous Runge-Kutta method with defect control, Numer. Alg., 16:349-364, 1997.

• Enright, W.H. and Hayashi, H., Convergence analysis of the solution of retarded and neutral differential equations by continuous methods, SIAM J. Numer. Anal., 35:572-585, 1998.

• Ermentrout, B., Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Researchers and Students, SIAM, Philadelphia, USA, 2002.

• Hoppensteadt, F.C. and Jackiewicz, Z. , Numerical solution of a problem in the theory of epidemics, Appl. Numer. Math., 56:533-543, 2006.

• Jackiewicz, Z. and Lo E., Numerical solution of neutral functional differential equations by Adams methods in divided difference form, J. Comput. Appl. Math., 18:592-605, 2006.

• Martin, A. and Ruan, S., Predator-prey models with delay and prey harvesting, J. Math. Biol., 43:247-267, 2001.

• Neves, K.W., Automatic integration of functional differential equations: an approach, ACM Trans. Math. Softw., 1:357-368, 1975.

• Neves, K.W. and Feldstein A., Characterization of jump discontinuities for state dependent delay differential equations, J. Math. Anal. and Appl., 5:689-707, 1976.

• Neves, K.W. and Thompson S., Software for the numerical solution of systems of functional differential equations with state dependent delays, Appl. Numer. Math., 9:385-401, 1992.

• Paul, C.A.H., A user-guide to ARCHI, Numerical Analysis Report 283, Mathematics Department, University of Manchester, U.K., 1995.

• Raghothama, A. and Narayanan, S., Periodic response and chaos in nonlinear systems with parametric excitation and time delay, Nonlin. Dyn., 27:341-365, 2002.

• Shampine, L.F., Numerical Solution of Ordinary Differential Equations, Chapman & Hall, New York, NY, 1994.

• Shampine, L.F., Solving ODEs and DDEs with residual control, Appl. Numer. Math., 52:113-127, 2005.

• Shampine, L.F., Gladwell, I., and Thompson, S., Solving ODEs with MATLAB, Cambridge Univ. Press, Cambridge, 2003.

• Shampine, L.F. and Thompson, S., Solving DDEs in Matlab, Appl. Numer. Math., 37:441-458, 2001.

• Thompson, S. and Shampine, L.F., A friendly fortran 90 DDE solver, Appl. Numer. Math., 56:503-516, 2006.

• Willé, D.R. and Baker, C.T.H., DELSOL - a numerical code for the solution of systems of delay-differential equations, Appl. Numer. Math., 9:223-234, 1992.

Internal references

• Zdzislaw Jackiewicz (2007) General linear methods. Scholarpedia, 2(4):2852.
• Frank Hoppensteadt (2006) Predator-prey model. Scholarpedia, 1(10):1563.
• John Butcher (2007) Runge-Kutta methods. Scholarpedia, 2(9):3147.
• Bard Ermentrout (2007) XPPAUT. Scholarpedia, 2(1):1399.

External Links

• Skip Thompson's website
• archi: http://www.maths.manchester.ac.uk/~chris/index.html
• dde-biftool: http://www.cs.kuleuven.ac.be/cwis/research/twr/research/software/delay/ddebiftool.shtml
• dde_solver: http://www.radford.edu/~thompson/ffddes/index.html
• dimsim: http://math.la.asu.edu/~jackiewi/zdzislaw.html
• dklag6: http://www.radford.edu/~thompson/misc/dklag.html
• radar5: http://www.unige.ch/~hairer/software.html
• snddelm: http://www.susqu.edu/facstaff/l/lo/

See Also

Numerical Methods, Ordinary Differential Equations

Sponsored by: Eugene M. Izhikevich, Editor-in-Chief of Scholarpedia, the peer-reviewed open-access encyclopedia

Reviewed by: Anonymous

Reviewed by: Dr. Zdzislaw Jackiewicz, Department of Mathematics, Arizona State University, Tempe, AZ

Accepted on: 2007-03-05 03:11:03 GMT