avoiding neg values in NDSolve
- To: mathgroup at smc.vnet.net
- Subject: [mg27317] avoiding neg values in NDSolve
- From: Michael Gilchrist <mag5 at duke.edu>
- Date: Tue, 20 Feb 2001 03:05:27 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
Hi all,
I've got a numerical issue that I would appreciate some advice on.
I am trying to solve an array of ODE's that describe a stochastic birth
death process where the mortality rate increases over time. I run into
a problem where NDSolve starts returning negative values for the
probability of being in state P[i] at time t.
I pretty much grasp that these negative values come from the fact that
0.0 != 0, but I'm wondering if anyone has a good suggestion on how to
minimize this type of error. Given that I will at some point
have hundreds of these ODEs to solve, increasing WorkingPrecision seems
like it will slow everything to a halt.
I'm wondering if there is a way to use Chop or something like that.
I've included a copy of the routine I've written and any suggestions
would be greatly appreciated.
Thanks,
Mike
(*routine*)
SolveEqn[m_, r_, p0_, tstart_, tstop_, popmax_:100] :=
Module[{vec, ivec, P, \[Mu], t, i, eqnvec, ndsolvevec, tee, sol},
(*make list of initial distribution at tstart*)
vec = ReplacePart[Table[0, {i, 0, popmax}], 1, p0 + 1];
ivec = Table[P[i][0] == vec[[i + 1]], {i, 0, popmax}];
(*make list of the diffeqns*)
diffeqvec =
Join[ {P[0]'[t] == \[Mu][t] P[1][t]},
Table[P[i]'[t] ==
r(i - 1) P[i - 1][t] - (r + \[Mu][t])i P[i][t] + \[Mu][
t] (i + 1)P[i + 1][t], {i, 1,
popmax - 1}], {P[popmax]'[t] == - \[Mu][t] popmax P[10][t]}] ;
(*list of variables to be solved by NDSolve*)
eqnvec = Table[P[i][t], {i, 0, popmax}];
(* list of eqns inc. IC, for NDSolve to solve*)
ndsolvevec = Join[diffeqvec, ivec];
tee = tstop - tstart;(*\[CapitalDelta]t*)
sol = NDSolve[ndsolvevec /. {\[Mu][t] -> m (t + tstart)},
eqnvec,
{t, 0, tee}, MaxStepSize -> 0.0001, MaxSteps -> 10000,
WorkingPrecision -> 30];
Evaluate[eqnvec /. sol] /. {t -> tee}
]
(*Here's it being used*)
In[120]:=
SolveEqn[0.01, 0.79, 4, 0, 0.25, 20]
Out[120]=
\!\({{5.65402894395491`*^-15, 6.767883228316017`*^-11,
3.038059606954285`*^-7, 0.0006061785202291144`,
0.45374239373450787`,
0.32510483856575756`, 0.14562015602314116`, 0.05218378013548495`,
0.01636319223414593`, 0.004691240128001236`,
0.0012609026068701344`,
0.0003227661485404878`, 0.00007952366621119406`,
0.00001899947990896227`, 4.4258050280020575`*^-6,
1.0093349551336155`*^-6, 2.2607543554495168`*^-7,
4.985839904782445`*^-8,
1.084646958106019`*^-8, \(-8.64633012287939`*^-12\), \
\(-2.31638523016838`*^-6\)}}\)
Michael Gilchrist
Dept. of Biology
Duke University
Durham, NC 27708