如何解决R语言deSolve包ode命令中的DLSODA精度警告？

问题现象 #

运行ode代码时触发DLSODA精度警告，时间t=1时计算结果全部为NA：

• 错误提示：

  DLSODA- 在T(=R1)处，请求的精度超出机器精度限制，请查看TOLSF(=R2)
  上述信息中，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
  

• 警告信息：

  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

警告原因 #

1. 模型函数参数顺序错误：model1的定义顺序为(parameters,state,time)，但deSolve的ode函数要求模型函数的签名必须是(time, state, parameters)，顺序混乱会导致求解器传入参数错误，直接引发数值求解异常。
2. 参数与变量数量级差异过大：模型中参数跨度极大，比如gamma=1.4e-6、g=-0.33，状态变量B初始值为1e6，这种悬殊的数量级会让求解器在精度控制时超出机器处理范围。
3. 默认精度设置过严：ode默认的rtol=1e-6、atol=1e-6精度要求，在系统存在极端数量级差异时，超出了DLSODA求解器的处理能力，触发警告并提前终止计算。

解决办法 #

1. 修正模型函数参数顺序：按照ode要求的标准签名调整函数参数顺序为(time, state, parameters)。
2. 放宽精度限制：增大rtol和atol的值，比如设置为1e-4，适配系统的数量级差异。
3. 可选优化：变量归一化：对数量级极大的变量（如B）做归一化处理，比如除以1e6，计算完成后再还原，降低数值求解的压力。

修正后的代码示例 #

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)

# 修正参数顺序为time, state, parameters
model1 <- function(time, state, parameters){
  with(as.list(c(parameters, state)), {
    N_total <- S + I + R  
    dS <- (mu * N_total) - (beta * S * I / N_total) - 
      (gamma * B * S / (B + K)) - (mu * S) + (omega * R)
    
    dI <- (beta * S * I / N_total) + 
      (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 <- as.data.frame(ode(y = state, times = time, func = model1, parms = parameters, 
                            rtol = 1e-4, atol = 1e-4))
print(output)

内容的提问来源于stack exchange，提问作者PAC-MAN