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

```

• Prev by Date: Re: Capturing error messages
• Next by Date: Re: InitializationCellEvaluation
• Previous by thread: Re: Capturing error messages
• Next by thread: Re: avoiding neg values in NDSolve