首页 > 解决方案 > 如何调节 ODE 函数中的参数值?(deSolve, if-else)

问题描述

在给定状态的初始值的情况下,我正在尝试创建以某些参数为条件的值。例如,如果 D 状态为D >= 60,则 S 值将为 S=1800。否则,如果 D 状态为D <60,则 S 值为S=4800。我在 ode 函数 (AedesIbag_model) 中使用了函数if-else 。当我使用D=70运行 ode 时,if-else 不会切换 S 参数值。所以我还没有做到这一点。如果我的英语不是很好,我很抱歉。感谢您的任何帮助。

library(deSolve)

AedesIbag_model<-function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dL = R*theta*S - mu_L*L - gamma_L*L - mu_w*L  
    dP = gamma_L*L - mu_P*P  - gamma_P*P - mu_w*P   
    dA = gamma_P*P - mu_A*A 
    dD = beta - alpha*D 
    if (D >= 60) { 
      S = 1800
      } else if (D < 60) {
        S = 4800
        } else if (D >= 10) {
          mu_w = 0.1
          } else if (D < 60) {
            mu_w = 0.1*100
          }
    
    return(list(c(dL, dP, dA, dD)))
  })
}

parameters  <- list(R = 0.24, theta = 10, S = 0,
                gamma_L = 0.2088, gamma_P = 0.384,
                beta = 10, mu_L = 0.0105, mu_P = 0.01041,
                mu_A = 0.091, mu_w = 0.1, alpha = 10
                )
state      <- c(L = 100, P = 60, A = 10, D = 70)
times       <- seq(0, 100, by = 0.1)

out_1 <- ode(y = state, times = times, func = AedesIbag_model, parms = parameters)
parameters

当我运行我的模型时。条件参数不会改变值。看!!!

> parameters
$R
[1] 0.24

$theta
[1] 10

$S
[1] 0   #S value doesn't change

$gamma_L
[1] 0.2088

$gamma_P
[1] 0.384

$beta
[1] 10

$mu_L
[1] 0.0105

$mu_P
[1] 0.01041

$mu_A
[1] 0.091

$mu_w
[1] 0        #S value doesn't change

$alpha
[1] 10

如果我的英语不是很好,我很抱歉。感谢您的任何帮助。

标签: rif-statementode

解决方案


我猜你正在使用 pkg deSolve。查看您的代码,首先要注意的是,尝试通过打印 的值来评估函数的行为parameters是没有效率的。该ode函数仅接受这些值,但不修改它们,因为 R 是一种函数式语言,因此不应修改其参数。下一个要纠正的错误是 if-else 结构。它永远不会修改mu_w其当前形式的参数,因为任何一个都D >= 60可能D < 60S修改为 never mu_w

我不知道 if(.){.}else{.} 构造是否会起作用,并开始相信你的话,它失败了,直到我意识到我上面写的内容。所以我使用了一种不同形式的逻辑运算来代替级联if-else用布尔数学构造:

    S    <- (D >= 60)* 1800 + (D < 60) * 4800
    mu_w <- (D >= 10) * 0.1 + (D < 60)* 0.1*100

我还认为函数参数的名称ode应该与目标函数中的名称相匹配,所以应该是parms而不是parameters,但是当我进行更改时我没有看到行为上有太大的不同,所以也许参数是按位置传递的,而不是比名字。如果您想查看ode调用结果如何演变,绘制结果更有效:

 png(); matplot(out_1, pch=1:3)
       legend("topright", 4, unlist(dimnames(out_1)[2]),
              pch = 1:5, col = 1:5)
 dev.off()

关于调试的进一步说明:由于 的值parameters未更改,因此您需要在定义的函数中放置一个printorcat语句来监视命名参数的本地环境值的变化。

在此处输入图像描述

使用packageDescription('deSolve')一个了解到有一个 gitHub 托管网页,其中包含教程。在该页面上进一步了解到有一个关于该主题的帮助页面?forcings。我强烈推荐非常易读的文本“在 R 中求解微分方程”。还有一个动态模型邮件列表的链接:邮件列表:https ://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

在简要查看教程和“强制”的材料后,我怀疑包作者建议使用分段线性函数而不是具有状态变化的参数的不连续函数。


推荐阅读