[R] Using odesolve to produce non-negative solutions

[dave fournier]

If you didn't get this solved.

I have done parameter estimation with models
defined by ODE's where negative solutions are a problem
and one can only avoid them with great difficulty if the
standard explicit methods for solving the ODE are used. I found that
using implicit methods could be a great help.

For example in the exponential case

     dN/dt = -k*N

the simple finite difference approximation is

     N_{t+1}-N+t
     --------------  = -k*N_t ,  k>=0
        h

   or
     N_{t+1} = N_t -k*h*N_t

   and if  k*h gets too large N_{t+1} goes negative and you are in trouble.

   Consider instead the implicit formulation where the second
   N_t on the RHS is replaced by N_{t+1}  and one gets

     N_{t+1} = N_t/(1+k*h)

    which is correct  for k*h=0 and as k*h--> infinity

   For a more complicated example see

    http://otter-rsch.com/admodel/cc4.html

   for something I called "semi-implicit".
   I hope these ideas will be useful for your problem.


         Cheers,

          Dave





-- 
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-rsch.com