[R-sig-dyn-mod] System with an event leading to a discontinuity in the ODE equations: fear of non convergence and

[Setzer.Woodrow at epamail.epa.gov]

The problem is that the ode solvers assume that the state variables (S,
V, D, E, and I) are continuous and smooth with smooth derivatives.  Your
'ifelse' statement guarantees that they are not, since the derivatives
jump brom 0 to nu.  The way to do this is to break the integration into
parts where the smoothness conditions are met - at t == tres and S == 0.
Take a look at lsodar, and the new 'events' feature in the version of
deSolve that has just been uploaded to CRAN.  That will allow you to
define an event (a function that changes state variables in discrete
jumps) when one of a number of conditions are met.  The solver will
automatically stop and restart the integration at that point.

R. Woodrow Setzer, Ph. D.
National Center for Computational Toxicology
http://www.epa.gov/comptox
US Environmental Protection Agency
Mail Drop B205-01/US EPA/RTP, NC 27711
Ph: (919) 541-0128    Fax: (919) 541-1194


                                                                                                                  
  From:       Raffaello Vardavas <r_vardavas at hotmail.com>                                                         
                                                                                                                  
  To:         <r-sig-dynamic-models at r-project.org>                                                                
                                                                                                                  
  Date:       12/18/2009 12:20 PM                                                                                 
                                                                                                                  
  Subject:    [R-sig-dyn-mod] System with an event leading to a discontinuity in the ODE equations: fear of non   
              convergence and                                                                                     
                                                                                                                  
  Sent by:    r-sig-dynamic-models-bounces at r-project.org                                                          
                                                                                                                  






Dear
Dr. Soetaert,

My
name is Raffaele Vardavas and I work mainly on
infectious disease modeling.
I have
come across your excellent online source on the R package DeSolve. I
usually use
Mathematica for solving
my ODE
models but since I am getting more familiar with R, I would like to
transition
to using DeSolve package in R.

I have
however a small problem and I was wondering if you may kindly sparse
some time
to help me with.

I have
written some code that integrates the following model for smallpox
vaccination
dynamics specified by Bauch et al. in PNAS 2003

 SEIDV_model <- function(t, Omega, parameters)
{
 with(as.list(c(Omega,parameters)), {
   dV = ifelse(t<tres | S<=0
,0, nu)
   dS = -beta*S*I-dV
   dE =  beta*S*I-sigma*E
   dI =
sigma*E-gamma*I
   dD =  gamma*I
   return(list(c(dS, dE, dI, dD,
dV)))
   })
   }

As you
can see this set of ODEs have a condition: the vaccination rate dV/dt is
nonzero
as long as the susceptible population is positive and
vaccines have been made available (after time tres).

When I
integrate this model I get many warnings:

Warning messages:
1: In lsoda(y, times, func, parms, ...) :

DLSODA-  At current T (=R1), MXSTEP (=I1) steps


2: In lsoda(y, times, func, parms, ...) :
        taken on this call
before reaching TOUT
3: In lsoda(y, times,
func, parms, ...) :
        In above message,  I1 =
5000
4: In lsoda(y, times, func,
parms, ...) :
        In above message,  R1 =
0.2381220791874D+02
5: In lsoda(y, times, func,
parms, ...) :
  an excessive amount of work (> maxsteps ) was done, but
integration was successful - increase maxsteps
6: In lsoda(y, times, func,
parms, ...) :

I have
increased the option maxstep and played around with hmax and atol
options - but
I seem not to be able to get rid of these warnings.
I do
get results in agreement with the authors of the paper. However,
I am going to be making this model more complicated later on with
additional boundaries ( thus more if statements
in the
ODE equations )- so I would rather be
comfortable in knowing and managing these warnings before
proceeding.

Could
you kindly shed some light on this.


Your
help would be very much appreciated.

Thank
you.

Raff.


_________________________________________________________________
Keep your friends updated—even when you’re not signed in.

cial-network-basics.aspx?ocid=PID23461::T:WLMTAGL:ON:WL:en-xm:SI_SB_5:092

             [[alternative HTML version deleted]]

_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models