Re: avoiding neg values in NDSolve
- To: mathgroup at smc.vnet.net
- Subject: [mg27334] Re: avoiding neg values in NDSolve
- From: "Allan Hayes" <hay at haystack.demon.co.uk>
- Date: Wed, 21 Feb 2001 03:17:11 -0500 (EST)
- References: <96t9gv$prt@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Michael,
If you are only concerned with the output then this can be adjusted.
Here are two suggestions
Call your output sol:
sol={{5.65402894395492`*^-15,6.767883228316035`*^-11,3.0380596069542927`*^-7
,
0.0006061785202291156`,0.4537423937345088`,0.32510483856575734`,
0.14562015602314077`,0.05218378013548467`,0.01636319223414579`,
0.004691240128001182`,0.0012609026068701157`,0.00032276614854048167`,
0.00007952366621119239`,0.000018999479908961843`,
4.425805028001938`*^-6,1.0093349551335863`*^-6,
2.2607543554494464`*^-7,4.98583990478228`*^-8,
1.0846469581059816`*^-8,-8.646330122923448`*^-12,-2.316385230168342`*^\
-6}};
First way
Chop[sol ,10^-5]
{{0, 0, 0, 0.0006061785202291156, 0.4537423937345088,
0.32510483856575734, 0.14562015602314077,
0.05218378013548467, 0.01636319223414579,
0.004691240128001182, 0.0012609026068701157,
0.00032276614854048167, 0.00007952366621119239,
0.000018999479908961843, 0, 0, 0, 0, 0, 0, 0}}
Second way
sol/.x_?Negative->0.
{{5.65402894395492*^-15, 6.767883228316035*^-11,
3.0380596069542927*^-7, 0.0006061785202291156,
0.4537423937345088, 0.32510483856575734,
0.14562015602314077, 0.05218378013548467,
0.01636319223414579, 0.004691240128001182,
0.0012609026068701157, 0.00032276614854048167,
0.00007952366621119239, 0.000018999479908961843,
4.425805028001938*^-6, 1.0093349551335863*^-6,
2.2607543554494464*^-7, 4.98583990478228*^-8,
1.0846469581059816*^-8, 0., 0.}}
--
Allan
---------------------
Allan Hayes
Mathematica Training and Consulting
Leicester UK
www.haystack.demon.co.uk
hay at haystack.demon.co.uk
Voice: +44 (0)116 271 4198
Fax: +44 (0)870 164 0565
"Michael Gilchrist" <mag5 at duke.edu> wrote in message
news:96t9gv$prt at smc.vnet.net...
> 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
>
>
>
>