Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Delay Differential Equations

• To: mathgroup at smc.vnet.net
• Subject: [mg8890] Re: Delay Differential Equations
• From: sherod at newton.Colorado.EDU (Scott Herod)
• Date: Tue, 30 Sep 1997 20:17:02 -0400
• Organization: /usr/local/lib/rn/organization
• Sender: owner-wri-mathgroup at wolfram.com

In article <60e6lb$c24 at smc.vnet.net> Peter Waterman <waterman at math.niu.edu> writes:
>Can I solve DDE of the form
>
>x'(t) = F(t, x(t), x(t-a))
>
>with NDSolve? It appears it is not able to refer back to its own
>solution while it is still in progress.
>
>Matthias Weiss
>mweiss at math.niu.edu
>

I wrote a package a year or so ago that is supposed to mimic NDSolve but
work on delay-differential equations.  Basically it loops over a time
interval of the delay time and looks up the solution from the
previous loop.  (I also fixed it to allow multiple delay times.)

Below is the code:

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Options[NDelaySolve] = {Verbose -> False};

NDelaySolve[eqnConst_, initsConst_, indinterval_, opts___] :=
  Module[{
    error = "Something is wrong",
    ind = indinterval[[1]],
    start = indinterval[[2]],
    stop = indinterval[[3]],
    eqn = Flatten[{eqnConst}],
    inits = Flatten[{initsConst}],
    tau, deps, numdeps, dtymes, delayedfs, numb,
    numbdelay, delaynum, onedelay, forcedterms, sol,
    functs, witchlist, g, numiters, subs, senteqn,
    i,j,k1,k2},
  
  deps = Map[Part[#,1] &, inits];
  numdeps = Length[deps];
  
  printQ = Verbose /. {opts} /. Options[NDelaySolve];
  If[printQ,Print["eqn = ", eqn]];
  If[printQ,Print["inits = ", inits]];
  If[printQ, Print["interval = " , indinterval]];
  If[printQ, Print["ind = " , ind]];
  If[printQ, Print["start = " , start]];
  If[printQ, Print["stop = " , stop]];
  If[printQ, Print["deps = ", deps]];
  
  dtymes = Map[ (ind-#)&,
    Union[Flatten[delayedfs = Map[Cases[#, _Plus, Infinity]&,
      Map[Cases[eqn,#, Infinity] &, deps /. ind -> Plus[__]]]
    ]]];
  
  If[printQ,Print["dtymes = ", dtymes]];

  dtymes = Sort[
    Map[If[MemberQ[{Integer, Rational, Real}, Head[#]],#,
    If[Head[(numb = N[#])] === Real, numb, error]] &,dtymes]];
    
  If[printQ,Print["dtymes = ", dtymes]];

  If[MemberQ[dtymes,error],Return[error]];
  If[(tau = dtymes[[1]]) <= 0, Return[error]];
  onedelay = ((numbdelay = Length[dtymes]) == 1);

  delayedfs = Map[Union,delayedfs];
    If[printQ,Print["delayedfs = ", delayedfs]];
  
  forcedterms = Table[Table[deps[[i]] /. ind -> delayedfs[[i,j]],
    {j,1,Length[delayedfs[[i]]]}], {i,1,Length[delayedfs]}];
  If[printQ, Print["forcedterms = ", forcedterms]];
  
  delayedfs = N[delayedfs];
  delaynum = Map[Length, delayedfs];
  
  
  funcs = Map[Apply[Rule,#]&, inits];
  If[printQ, Print["funcs = ", funcs]];
  witchlist = Table[Which[Evaluate[ind <= start,
    deps[[i]] /. funcs[[i]]] ], {i,1,numdeps}];
  g = Map[Function[Evaluate[ind], #] &, witchlist];
  numiters = Ceiling[(stop - start)/tau];

  For[j = 1, j <= numiters,j++,
    jt = N[start + j*tau];
    If[printQ, Print["jt = ", jt]];
  
    sub = Flatten[Table[
      Table[forcedterms[[k2,k1]] -> g[[k2]][delayedfs[[k2,k1]]],
      {k1,1,delaynum[[k2]]}], {k2,1,Length[delaynum]}]];
   
    ics = Table[(deps[[i]] /. ind -> jt - tau) ==
      g[[i]][jt - tau], {i,1,Length[deps]}];
      
    senteqn = Union[eqn /. sub, ics];
    
    If[printQ, Print["senteqn ", j, " = ", senteqn]];

    sol = Flatten[
      NDSolve[senteqn, deps, {ind, jt - tau, jt},opts]];

    If[onedelay,
      witchlist = Table[Join[witchlist[[i]], 
        Apply[Which,{ind <= jt,
        deps[[i]] /. sol}]], {i,1,numdeps}];
      g = Map[Function[Evaluate[ind], #] &,
        Table[Apply[Which,{ind <= jt,
        deps[[i]] /. sol}], {i,1,numdeps}]],
      g = Map[Function[Evaluate[ind], #] &,
        Table[Join[g[[i,2]], Apply[Which,{ind <= jt,
        deps[[i]] /. sol}]], {i,1,numdeps}]]
    ]
  ];
  If[onedelay,
    g = Map[Function[Evaluate[ind], #] &, witchlist]];
        
{Thread[Rule[deps,Through[g[ind]]]]}

];

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Call it with something like:

sol = NDelaySolve[{x'[t] == - x[t] y[t-1] + y[t-10],
             y'[t] == x[t] y[t-1] - y[t],
             z'[t] == y[t] - y[t-10]},
  {x[t] == 5.0, y[t] == 0.1, z[t] == 1.0}, {t,0,100.0}];


Notice that this mimics the syntax of NDSolve but initial values
(which may be functions of t) are included in the dependent variable
specification section (the second argument).

Also note that the above example has two delay times, t0 = 1 and
t1 = 10.

Scott A. Herod
Applied Mathematics
University of Colorado, Boulder
Scott.Herod at Colorado.EDU