MathGroup Archive 2001

[Date Index] [Thread Index] [Author Index]

Search the Archive

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




  • Prev by Date: Re: chaos-to -order transform
  • Next by Date: Re: 2D Outlines to 3D Surface or Solid
  • Previous by thread: avoiding neg values in NDSolve
  • Next by thread: Re:avoiding neg values in NDSolve