Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2001
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2001

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

Search the Archive

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