How to tackle the warning in ode command of desolve package in R?

Question

I have tried to write this which is otherwise fine but I do not know what is the reason I am getting the error.

DLSODA-  At T (=R1), too much accuracy requested  
      for precision of machine..  See TOLSF (=R2) 
In above message, R1 = 1, R2 = nan
 
  time    S  I  R     B
1    0 9999  1  0 1e+06
2    1   NA NA NA    NA
3    1   NA NA NA    NA

This is the output I am getting.

I have tried to run this code and it doesnot look to have any problem for me .I am adding my code in the next line. Can anyone tell me the reason of this error ?

library(deSolve)  
t=30
beta=0.058
gamma=1.4*exp(-6)
N=10000
state=c(S=N-1,I=1,R=0,B=10^6)
time<-seq(from=0, to =t,by=1)
parameters <- c(mu=5.089041e-05,phi=0.2,g=-0.33,xi=10,K=1e6,omega=0.00137,
                beta=beta,gamma=gamma)
model1<- function(parameters,state,time){
  with(as.list(c(parameters,state)), {
    N=state['S']+state['I']+state['R']  
    dS=(parameters['mu']*N)- (parameters['beta']*state['S']*state['I']/N)- 
      (parameters['gamma']*state['B']*state['S']/(state['B']+parameters['K']))-(parameters['mu']*state['S']) +
      (parameters['omega']*state['R'])
    
    dI=(parameters['beta']*state['S']*state['I']/N) + 
      (parameters['gamma']*state['B']*state['S']/(state['B']+parameters['K']))-(parameters['mu']*state['I'])   -
      (parameters['phi']*state['I'])
    
    dR=(parameters['phi']*state['I'])-(parameters['mu']*state['R']) -(parameters['omega']*state['R'])
    
    dB=(parameters['xi']*state['I'])-(parameters['g']*state['B'])
       
    return(list(c(dS, dI, dR, dB))) 
  })
}

output <- as.data.frame(ode(y = state,  times = time,func = model1, parms = parameters ))

The warning messages are as follows :

Warning messages:
1: In lsoda(y, times, func, parms, ...) :
  Excessive precision requested.  scale up `rtol' and `atol' e.g by the factor 10
2: In lsoda(y, times, func, parms, ...) :
  Returning early. Results are accurate, as far as they go

Any kind of help would be grateful. Thanks in advance.

Answer 1

The main error was, that the order of arguments in the model1-function was wrong, it needs to be always in a fixed order time, state, parms as written in the help file. This fixed order is necessary for performance reasons.

Besides this, I simplified the style to improve readability. If the with()-function is used anyway, names of the parameter and state variables can be directly accessed, as long as there are no duplicates. In addition, I replaced the assignment = to <- and removed the as.data.frame in the interest of the built-in plot method.

Please check yourself, whether gamma should be 1.4*exp(-6) (0.00347) or 1.4e-6 (0.0000014).

library(deSolve)  
t     <- 30
beta  <- 0.058
#gamma <- 1.4*exp(-6)
gamma <- 1.4e-6
N     <- 10000

state=c(S=N-1, I=1, R=0, B=1e6)

time<-seq(from=0, to = 30, by=1)

parameters <- c(mu=5.089041e-05,phi=0.2,
                g=-0.33,xi=10,K=1e6,omega=0.00137,
                beta=beta,gamma=gamma)

model1 <- function(time, state, parameters){
  with(as.list(c(parameters, state)), {
    
    N <- S+I+R  
    dS <- mu*N - beta*S*I/N - gamma*B*S/(B+K) - mu*S + omega*R
    dI <- beta*S*I/N + gamma*B*S/(B+K) - mu*I - phi*I
    dR <- phi*I - mu*R - omega*R
    dB <- xi*I - g*B
    
    return(list(c(dS, dI, dR, dB))) 
  })
}

output <- ode(y = state,  times = time,func = model1, parms = parameters)

plot(output)