Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
1997
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 1997

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

Search the Archive

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


  • Prev by Date: Re: Undefined macro in LaTex file from Mma 3.0
  • Next by Date: Export of Mathematica surfs to 3D file formats
  • Previous by thread: Re: Undefined macro in LaTex file from Mma 3.0
  • Next by thread: Export of Mathematica surfs to 3D file formats