2

我有一个微分方程模型（见下文），我正在寻找某种方法来改变某个时间步长或某个状态变量的某个值的某些参数值（或另一种方法，如果它工作得更好）。例如，我想在时间步 5 将 GammaQR 和 GammaQD 从 0 更改为 .02（或者，如果效果更好，当 H > .04 时可能）。我不知道该怎么做，非常感谢任何建议！谢谢！

 VS  = 0.01    
VE  = 0.01  
VH  = 0.01  
BSE = 10 #
BSI = 10  #
THE = 1/10      
TEI =  1/3 
DeltaID = .5 
GammaIR = .5 
XIQ = 0.001  
XEQ = 0.01 
XHQ = 0.05 
GammaQR = .02  
GammaQD = .02 

library(deSolve)

pars <- c(VS, VE, VH, BSE, BSI, THE, TEI, DeltaID, GammaIR, XIQ, XEQ, XHQ, GammaQR)

init.values <- c(S = .99 , H = .01 , E = 0 , V = 0, I = 0 , Q = 0 , D = 0 , R = 0 )

times <- seq(0, 60, by = 1)

Smallpox <- function(time, y.values, parameters){
   with(as.list(c(y.values, parameters)), {

    dS.dt = -VS*S - BSI*S*I - BSE*E*S 
    dH.dt = BSI*S*I + BSE*E*S  - THE*H - XHQ*H - VH*H 
    dE.dt = THE*H - XEQ*E - VE*E - DeltaID*E - TEI*E
    dV.dt = VS*S + VE*E + VH*H
    dI.dt = TEI*E - XIQ*I - GammaIR*I - DeltaID*I
    dQ.dt = XHQ*H + XIQ*I + XEQ*E - GammaQR*Q - GammaQD*Q
    dD.dt = DeltaID*I + DeltaID*E + GammaQD*Q
    dR.dt = GammaQR*Q + DeltaID*I


    return(list(c(dS.dt, dH.dt, dE.dt, dV.dt, dI.dt, dQ.dt, dD.dt, dR.dt)))
  })
}

out <- as.data.frame(ode(func = Smallpox, y = init.values, parms = pars, times = times))
tail(out)

matplot(out$time, out[ ,2:9], type = "l", xlab = "time", ylab = "percent of       population", main = "Model Name", lwd = 2, col = c("black", "red", "green", "Blue", "cyan", "purple", "grey", "magenta"), lty = 1:8)
4

0 回答 0