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

MathGroup Archive 2004

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

Search the Archive

Re: System of DE's

  • To: mathgroup at smc.vnet.net
  • Subject: [mg48230] Re: System of DE's
  • From: sean_incali at yahoo.com (sean kim)
  • Date: Wed, 19 May 2004 02:42:05 -0400 (EDT)
  • References: <c83t0e$mar$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

if you have Scott herod's NDelaySolve code then below works.


<<NDelaySolve.m

vars = {x[t], y[t], z[t]};
 k1 = 0.0001; 
k2 = 0.1; 
k3 = 0.0001; 
n = 2000; 

f1 = -k1 x[t] y[t] - k3 x[t] z[t]; 
f2 = -k2 y[t] + k1 x[t] y[t] + k3 x[t] z[t]; 
f3 = k2 (y[t] - y[t - 14]); 

inits = {x[t] == n - 10, y[t] == 10, z[t] == 0}; 
eqns = {x'[t] == f1, 
    y'[t] == f2, z'[t] == f3}; 

sol = NDelaySolve[ eqns,  inits,  {t, 0, 200}]; 

Plot[Evaluate[{x[t], y[t], z[t]} /. sol], {t, 0, 200}, PlotRange ->
All];


below is his code. I made a package loadable in Mathematica 5.0 using his
codes.
copy from "BeginPackage" to "EndPackage" at the end then paste intot
anotebook and asave as a package format. then put it into C:\Program
Files\Wolfram Research\Mathematica\5.0\AddOns\Applications folder

then <<NDelaySolve.m can call it.

sean

---------------


BeginPackage[ "NDelaySolve`"]

NDelaySolve::usage = "sol = NDelaySolve[
{x'[t] == - x[t] y[t - 1] + y[t - 6], y'[t] == x[t] y[t - 1] - y[t],
z'[t] == y[t - 1] - y[t - 6]},
{x[t] == 20.0, y[t] == 0.1, z[t] == 1.0}, 
{t, 0, 20.0}]; \n
\n
Plot[Evaluate[{x[t], y[t], z[t]} /. sol], {t, 0, 20}, PlotRange ->
All];
\n
\n
NDelayDSolve[eqns, initcond, {t,tmin,tmax}]\n 
\n
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).
\n
Also note that the above example has two delay times, t0 = 1 and t1 =
10.
\n
Scott A. Herod\n
Applied Mathematics\n
University of Colorado, Boulder\n
Scott.Herod at Colorado.EDU " 

Begin["`Private`"]

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]]]]} ]; 




End[]
EndPackage[]


  • Prev by Date: Re: From C/C++ to Mathematica?
  • Next by Date: Re: From C/C++ to Mathematica?--> MathLink programs in C for Windows
  • Previous by thread: Re: System of DE's
  • Next by thread: Simple, syntactical question