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

[Soetaert, Karline]

Raffaello,
 
 
By definition, a parameter does not change, so in deSolve, parameters always have the same value. 
 
To allow "nu" to change value after the event, you can make it a STATE VARIABLE, whose rate of change = 0. With its derivative = 0,  it will not change, except during the event.
 
Something like this (where Omega now =c(S,E,I,D,V,Nu)):
 
 
SEIDV_model <- function(t, Omega, parameters) {
 with(as.list(c(Omega,parameters)), {
   dV =  Nu
   dS = -beta*S*I-dV
   dE =  beta*S*I-sigma*E
   dI =  sigma*E-gamma*I
   dD =  gamma*I
   dNu = 0.
   return(list(c(dS, dE, dI, dD, dV, dNu)))
   })
   }
   ...
   
   # sets vacination rate nu to 0.
eventfun <- function(t, Omega, parameters) {
  return(list(c((Omega[1:5], Nu=0)))


This should work.
 
 
Karline
 
________________________________

From: r-sig-dynamic-models-bounces at r-project.org on behalf of Raffaello Vardavas
Sent: Fri 1/8/2010 12:12 AM
To: r-sig-dynamic-models at r-project.org
Subject: Re: [R-sig-dyn-mod] System with an event leading to a discontinuity in the ODE equations: fear of non convergence and




Dear Dr. Setzer,

thank you for your reply. I am exploring the use of lsodar and the 'events' option. I see that this is indeed the way to proceed.
However I am a little confuse on how to implement this.

my first attempt is shown below:

#############################################################################
##          SEIDV Epidemic Model
#############################################################################

SEIDV_model <- function(t, Omega, parameters) {
 with(as.list(c(Omega,parameters)), {
   dV =  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)))
   })
   }
  
# event triggered if state variable S<=0 or V>V.lim
rootfun <- function (t,  Omega, parameters) {
with(as.list(c(Omega,parameters)),{
  ifelse(S<=0 | V>V.lim, r <- 0 , r <- 1 )
return(r)})}

# sets vacination rate nu to 0.
eventfun <- function(t, Omega, parameters) {
with(as.list(c(Omega,parameters)),{
return(nu=0)})}


this however gives an error when I integrate with lsodar:
Error in checkEventFunc(Eventfunc, times, y, rho) :
  The number of values returned by events$func() (1) must equal the length of the initial conditions vector (5)

I believe this is because the eventfun requires that I specify the values of the new set of state variables (all 5 of them)
once the event condition is met and not a change of parameters.
However, I do not want to change the state variables but rather one parameter (nu) at this event.
(1) How can I do this?

Alternatively - the other way of proceeding would be to specify a second function
SEIDV_model_2 with the first equation set to dV=0. Thus the eventfun would need to stop the integration at the event and output the last set of state variables - from which I can then
use as initial conditions to SEIDV_model_2.

(2) How do I specify in eventfun that the integration should stop here.


Thank you for your kind help in this.

Best.

Raffaello Vardavas.




                                         
_________________________________________________________________


cial-network-basics.aspx?ocid=PID23461::T:WLMTAGL:ON:WL:en-xm:SI_SB_1: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


-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/ms-tnef
Size: 7658 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-dynamic-models/attachments/20100108/5e71ac8f/attachment-0001.bin>