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

[Thomas Petzoldt]

Dear Raffaelo,

thanks for your interest. You are right,  eventfun requires that you 
*specify* the values of the new set of state variables (all 5 of them).

However, it is of course not necessary that you *change* all of them. 
Simply return the old states as they are, e.g.:

# sets vaccination to 0 (not rate nu !!!)
eventfun <- function(t, Omega, parameters) {
   with(as.list(c(Omega,parameters)),{
     return(list(c(dS, dE, dI, dD, dV=0))
   })
}

# or alternatively:
eventfun <- function(t, Omega, parameters) {
   return(list(c((Omega[1:4], dV=0)))
}


Note also, that the event specifies an "immediate" drop of dV at the 
time of the event, so this requires to be passed directly to the state 
(dV) and not to the vaccination rate "nu" (that has a unit of 1/time).

If you instead want to set only vaccination *rate* nu to zero and let V 
unchanged, than it's another topic (a forcing function).

Short question: am I right that state variable V is the number of 
vaccinated people and do you want to set it to zero, because of an 
assumed mutation of a virus?

Please note that I could not test the model because I don't know the 
background of your model and you did not provide a *fully* reproducible 
example.

Nevertheless, hope it helps

Thomas Petzoldt


Raffaello Vardavas wrote:
> 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
> 


-- 
Thomas Petzoldt
Technische Universitaet Dresden
Institut fuer Hydrobiologie        thomas.petzoldt at tu-dresden.de
01062 Dresden                      http://tu-dresden.de/hydrobiologie/
GERMANY