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

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

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[
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