From: spellucci@mathematik.th-darmstadt.de (Peter Spellucci) Newsgroups: sci.math.num-analysis,sci.physics.research,sci.math.research,sci.nonlinear Subject: Re: Numerical Methods to integrate time-delay ODEs ? Date: 26 Feb 1998 11:07:56 -0800 In article , pecora@zoltar.nrl.navy.mil (Louis M. Pecora) writes: |> I am looking for information on the best numerical algorithms to use in |> integrating time-delay ordinary differential equations, e.g. |> |> dx |> -- = f(x(t), x(t-tau)) |> dt |> |> where the vector field depends on the dynamical variables at an earlier |> time. Is it sufficient to simply discretize the time interval [t-tau,t] |> and integrate a finite (but maybe large) number of ODEs or are there |> better, more accurate, and more efficient approaches? |> |> Thanks for any pointers, references, book sections, etc. There are lots of good methods known for this task, see e.g. in Hairer&Norsett&Wanner Solving ordinary differential equations I,II (Springer) There is code "retard.f" (from their book) in the library "hwlib" of the library "codelib" of elib.zib-berlin.de (http://elib.zib-berlin.de). This code is an adaptation of the Dormand-Prince (4)5 Runge-Kutta scheme. Similarly is 497 in netlib/toms. hope this helps peter ============================================================================== [Pecora then submitted collected responses when a similar question emerged-djr] ============================================================================== From: pecora@zoltar.nrl.navy.mil (Louis M. Pecora) Newsgroups: sci.math.num-analysis,sci.physics.research,sci.math.research,sci.nonlinear Subject: Re: Numerical Methods to integrate time-delay ODEs ? Date: 5 Mar 1998 09:26:28 GMT In article , Post.To.NewsGroup.NOT@email.thanks (Troy Shinbrot) wrote: > I do not have good pointers other than the ones you already know -- e.g. > Mackey-Glass. The key problem as I understand it is that the equation > becomes infinite dimensional when a delay is introduced, and one needs to > define an infinite dimensional set of initial conditions between t = 0 and > t = -tau. Then unless the particular problem is globally stable to > variations of these IC's, the behavior will depend on the set of IC's. > There may be a fantastic reference out there, which I would be glad to > see, but my understanding is that this is a fundamental limitation to our > understanding of this class of prob's. ============================================================================== [I am raising all nested quotes to top level; ] [ I hope this preserves attributions correctly --- djr. ] Assuming you have in mind an initial-value problem (i.e. the values x(t) are known for t0 - tau <= t <= t0 and you want to solve for x(t) where t > t0), you can use one of the usual methods for integrating (non-time delayed) IVP's, e.g. Runge-Kutta or predictor-corrector methods, with minor modifications. You will have to save the values of your approximate solution (at a discrete set of times, of course) on the interval from t - tau to t (and update this as your algorithm advances t), and use this saved information in the evaluation of the RHS of the ODE. Using "stupid programming tricks" you should even be able to use a canned routine for solving non-delayed ODE's of the form dx/dt = F(t, x(t)): look for a routine which has the calling sequence ODE_Solve(t, h, x, F) and which given t, a time step h, x(t), and (a program to calculate) F(t, x), returns an approximation to x(t + h). (This sort of solver is pretty common.) You can use such a solver to compute an array of approximate values for x at discrete times; step to the value at the next time by calling the solver with the current time, desired step, current x, and a function F which gets hold of the array of computed past values for x (which you have to save, of course), either as an additional parameter (if your solver permits this; most should) or, if necessary, as "external static" data (in C) or in a common block (in Fortran), and which uses this array to find x(t - tau) and then returns f(x, x(t - tau)). This kind of approach should be much more effective than approximating with lots of non-delayed ODE's. If you need more details, let me know. Robert E. Beaudoin ============================================================================== Check http://www.unige.ch/math/folks/hairer/software.html ============================================================================== From: Christopher T H Baker Date: Thu, 21 Apr 94 16:21:56 BST Subject: Delay Differential Equations NA-NET recently carried a request for references on delay differential equations from Heru Suhartanto (heru@cs.ui.ac.id). This prompts us to publicise the availability of our Numerical Analysis Report 248, Issues in the Numerical Solution of Evolutionary Delay Differential Equations by C T H Baker, C A H Paul and D R Wille available by anonymous ftp from vtx.ma.man.ac.uk (130.88.16.2) in pub/narep as narep248.ps.Z. This report is based on the talk at SCADE93 in Auckland, given by Christopher Baker, and is dedicated to John Butcher. It includes 75 references to existing work in this area. Incidentally, the authors would appreciate receiving comments on this report from others working in the area. Christopher Baker ============================================================================== From: mroussel@alchemy.chem.utoronto.ca (Marc Roussel) Newsgroups: sci.nonlinear Subject: Re: Mackey-Glass equation Date: 6 Jul 93 00:20:58 GMT In article diek@cs.utexas.edu (Diek Winters Wheeler) writes: >I want to generate a time series based on the delay-differential Mackey-Glass >equation. Can someone explain to me how one goes about numerically integrating >this equation? I have a program that integrates systems, such as the Lorenz >equations, with a fourth-order Runga-Kutta algorithm. Will an RK4 algorithm >also work with a delay-differential equation? Yes, it will, but you have to arrange to store the solution as it's generated and you have to provide an interpolator (of some sort) to provide the integrator with delayed variables as they are needed. This is a little tricky, but certainly not impossible. You should also be very careful. DDEs are somewhat more prone to numerical instability than ODEs so you should perturb things like the step sizes between runs to see that the results are reproducible. (This is probably a good idea for ODE systems as well, but ODE integrators are getting so reliable that this often feels like a waste of time.) Just to familiarize yourself with some of the issues, you might want to read some of the following papers: George Marsaglia, Arif Zaman and John C.W. Marsaglia, Math. Comp. 53, 191-201 (1989). This paper describes a very elegant method of solution based on power series expansions between the knots. G.S. Virk, IEE Proc. D 132, 119-123 (1985). This paper extends Runge-Kutta methods to delay-differential equations. The method presented is probably a reasonable compromise between computational efficiency and code complexity. Irving R. Epstein and Yin Luo, J. Chem. Phys. 95, 244-254 (1991). In addition to the application of delay-differential equations to kinetics, several simple numerical integration methods are presented. Marc R. Roussel mroussel@alchemy.chem.utoronto.ca ============================================================================== Newsgroups: sci.math.num-analysis Subject: Re: Solving the Mackey-Glass Chaotic Time-Series Date: Fri, 02 Feb 96 12:44:10 MEZ Organization: RHRZ Uni-Bonn In article bushbo@bobcat.ent.ohiou.edu (Brian O. Bush) writes: >I am attempting to numerically solve the Mackey-Glass >Chaotic Time-Series: > >dx(t) 0.2*x*(t-tau) >_____ = _____________ - 0.1*x(t) > dt 1 + x^10(t-tau) > > >with tau = 30 (the t.s. gets chaotic with tau > 17), it has roots in >the article mentioned below. > >I have read in several sources that a fourth-order Runge-Kutta method >is used to solve this. However, since this is a delay equation, we know that >at t = 0, x(t) = 0.8, but we need x(t-tau), which is much less than zero. >Thus, there needs to be some "bootstrap" method to solve this in maybe a >history function or something. > >I would be interested in how i should approach this, with the RK method. I >have also found little references on delay equations and their solutions >(and approaches). Thanks. Delay equations have some properties that distinguish them from 'simple' ordinary differential equations. E.g. You need to know not only the start value x(0)=x0, but a start FUNCTION x(t)=x0(t) for -tau <= t <= 0. If this start function does not fit to the DE (and usually it doesn't) then there is a jump discontinuity at t=0 either of x or some derivative D^n x. At t=tau (where (t-tau) 'crosses' 0) again there is a jump (usually the jump occurs now in the (n+1)st derivative), a next jump at t=2*tau and so on, each time the jump occuring one derivative 'later'. Runge-Kutta-Methods are only appropriate for SMOOTH solutions. But since x is smooth between the jumps, it is o.k. to use RK-methods as long as you don't integrate 'across' a jump discontinuity, i.e. you have to include the critical points into the RK-mesh. In your case all critical points are known in advance (0, tau, 2*tau, 3*tau, ...) but in general the delay may be a function of t or even of x(t) and t. In this case you need elaborate methods for keeping track of the jump point propagation. If you try to solve the DODE (delayed ODE) numerically via a RK-method you need to evaluate the solution x not only at the current time t, but also at some delay time (t-tau). The basic RKM will provide x only at times t that have been included to the mesh (0, h1, h1+h2, h1+h2+h3, ..., where h1, h2, ... are the step sizes, usually automatically chosen by the RKM). In your case, the points tau, 2*tau, etc. are (artificially) added to this mesh, because of the jump points, but that does NOT mean that when your integrator chooses some stepsize h and with this a new time t, you have already computed x(t-tau). Indeed, this only happens in very exceptional cases, namely if h=tau is chosen. There are two possible solutions to this: 1. You interpolate x(t-tau) using the neighboring mesh points t1 < t-tau < t2 and the known values x(t1), x(t2), x'(t1), x'(t2) (the derivative is known via the ODE) and performing a Newton inter- polation. 2. You use a 'better' RKM, namely a RK-Sarafian method or a RK-Fehlberg method, that provide not only a numerical solution at the mesh points, but also a (polynomial) solution BETWEEN these points. Remark: Of course the stepsize h must never exceed tau! To evaluate x(t-tau), the last meshpoint must lie between t-tau and t. Maybe you wish to use some standard software to solve your problem instead of creating your own. I recommend you to read Neves & Thompson "Software for the numerical solution of systems of functional differential equations with state dependent delays", J. Applied Num. Math., Vol. 9, 1992, pp. 385-401. If you are interested in their FORTRAN code for the program DRKLAG (which is discussed in the quoted article), let me know. But I have to admit that I never used the program and only studied it theoretically. I hope this helps a bit. Arndt. ============================================================================== Hairer, Norsett and Wanner, Solving Ordinary Differential Equations, has a section on delay differential equations. I use this book (actually, a 2-vol set) as a reference for my course at Cornell and recommend it highly. -- Steve Vavasis ============================================================================== Lou, that is a very good question. For linear TD-ODE's one can use transform (Laplace, Fourier, ...) methods, and something is known. For the nonlinear case that we'd like to look at, certainly the discretization approach has been used a lot. I think it may miss interesting phenomena. One of the problems in such studies is that one must specify an initial function x(t)= c(t) on [-tau, 0]. Let's suppose that the initial function c() is polynomial in t, and that f() is polynomial; one might then look for Taylor series expansions. If c() is a trig polynomial, one might look for periodic solutions x() of period tau... or look for almost-periodic solutions. I seem to remember that 30 years ago the almost-periodic solutions were of interest to specialists. It may be that Jack Hale would help... the Editor of J. Diff. Eqs., now at Georgia Tech. (Hale@math.gatech.edu) ... he is older and wiser than me. The applications have been in population dynamics, so you should look at the time-delay versions of Volterra-Lotka systems and other ecological models. There must be hundreds of papers in that literature. David David L. Elliott, Vis. Sr. Research Scientist, Inst. for Systems Research, ==============================================================================