R中deSolve用approxfun时，ODE系统时变参数下累积指标异常问题

在带年度时变参数的确定性仓室模型中，统计存活总人数、死亡总人数等累积指标时，使用approxfun处理时变参数后，添加累积指标会改变原有指标的统计结果。

以含易感者（S）、感染者（I）的SI模型为例：I状态的死亡率m每年变化，通过approxfun获取对应时间点的死亡率。初始模型（无死亡累积状态M）的存活总人数为5960657；添加死亡累积状态M后，存活总人数变为5960676。参数非时变时无此问题，原模型中差异更显著。

初始模型代码（无死亡累积状态） #

library(deSolve)

dt <- 0.5
times <- seq(from = 0, to = 100, by = dt)

# 每年死亡率变化——用均匀分布随机值模拟
set.seed(657)
mort_year <- runif(n=max(times), min=0, max=0.1)
# 使用approxfun获取对应时间的死亡率
mort <- approxfun(mort_year, method = "constant", rule = 2)

# 各状态初始人数
S_0 = c(999)                    
I_0 = c(1)    

si_initial_state_values <- c(S = S_0,
                             I = I_0)

# 参数值
si_parameters <- c(beta = 0.08, b = 0.05) 

# 模型定义
si_model <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    S <- state[1]
    I <- state[2]
    
    # 存活总人口
    N <- S + I
    
    # 感染力
    lambda <- beta * I / N  
    
    # 死亡率
    m <- mort(time)
    
    # 微分方程求解
    dS <- b*N -lambda * S 
    dI <-  lambda * S - m*I
    
    list(c(dS, dI))
  })
}

output <- ode(y = si_initial_state_values,
              time = times,
              func = si_model,
              parms = si_parameters)

N <- output[,2:3] 
sum(N)
# [1] 5960657  

添加死亡累积状态后的代码 #

# 各状态初始人数
S_0 = c(999)                    
I_0 = c(1)    
M_0 = c(0)

si_initial_state_values <- c(S = S_0,
                             I = I_0,
                             M = M_0)

# 参数值
si_parameters <- c(beta = 0.08, b = 0.05) 

# 模型定义
si_model <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    S <- state[1]
    I <- state[2]
    M <- state[3]
    
    # 存活总人口
    N <- S + I
    
    # 感染力
    lambda <- beta * I / N  
    
    # 死亡率
    m <- mort(time)
    
    # 微分方程求解
    dS <- b*N -lambda * S 
    dI <-  lambda * S - m*I
    dM <- m*I
    
    list(c(dS, dI, dM))
  })
}

output <- ode(y = si_initial_state_values,
              time = times,
              func = si_model,
              parms = si_parameters)

N <- output[,2:3] 
sum(N)
# [1] 5960676 

内容的提问来源于stack exchange，提问作者Andyga