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