R中的deSolve包，变量未更新

Question

我正在使用R和deSolve软件包来解决一组微分方程。变量为WC和C(土层中的含水量和浓度)。土层数量可以在代码中更改，如下所示：

numboxes <- 1                   # Number of soil layers
delx <- rep(1,numboxes)     # Thickness of soil layers (cm)
delt <- 1   
bulk<-0.5;FC<-0.4;Ks<-0.03;Sat<-0.8;Wres<-0.1;kbio=0.01;kd=10  #parameters
## the model
SPEC <- function(t,state,parms) {
    with(as.list(c(state,parms)),{
    ifelse (WC>=Wres, perc <- (WC*delx*10*bulk-FC*delx*10*bulk)*(1-exp(
    -Ks/(Sat*delx*10*bulk-FC*delx*10*bulk))), perc <-0)
    #
    dWC <- -diff(c(0,perc))*24/(10*delx*bulk)
    dC <- -kbio*C
    list(c(dWC,dC),perc=perc)
    })
}

然后使用deSlove包求解该函数：

WC <-rep(0.5,times=numboxes)
C <- rep(20,times=numboxes)
state <- c(WC=WC,C=C)
times <- seq(from=1, to=5, by=delt)
out <- as.data.frame(ode(times=times,y=state,func=SPEC,parms=0,method="rk4"))

对于单个土层(numboxes=1)，一切工作正常，结果如下：

  time        WC        C        perc
1    1 0.5000000 20.00000 0.007444030
2    2 0.4699599 19.80100 0.005207836
3    3 0.4489439 19.60397 0.003643397
4    4 0.4342411 19.40891 0.002548917
5    5 0.4239550 19.21579 0.001783219

但是，当我增加土层的数量(例如numboxes=2)时，求解器运行，但结果不正确。对于两个层：

  time       WC1 WC2   C1   C2      perc1      perc2
1    1 0.5000000 0.5 20.0 20.0 0.00744403 0.00744403
2    2 0.4642687 0.5 19.8 19.8 0.00744403 0.00744403
3    3 0.4285373 0.5 19.6 19.6 0.00744403 0.00744403
4    4 0.3928060 0.5 19.4 19.4 0.00744403 0.00744403
5    5 0.3570746 0.5 19.2 19.2 0.00744403 0.00744403

两个土层(C1和C2)的浓度计算结果是正确的，并与单层浓度的结果一致。但是，计算的含水量不正确(第一层的结果应与使用单层的模拟结果相同)。它认为，为了计算渗流(perc)，求解器仅使用WC0.5的初始值，因此每次迭代都会计算相同的值(perc1和perc2始终等于0.007，这对于WC=0.5是正确的)。

奇怪的是，当我使用单层时，情况并非如此。这个问题看起来与这里报道的类似：Population values not updating in deSolve in R然而，我正在以"state=value“的方式定义初始值，并且我仍然面临更新问题。你知道我该如何解决这个问题吗?为什么我要面对这个问题？

Answer 1

尝试将赋值放在ifelse之外的perc中。

perc <- ifelse(WC >= Wres,
               (WC*delx*10*bulk - FC*delx*10*bulk)*
               (1 - exp(-Ks/(Sat*delx*10*bulk - FC*delx*10*bulk))),
               0)

我甚至不知道在ifelse中做赋值会做什么。

Answer 2

您的问题是，您错误地定义了状态变量。运行代码，然后

state

你看，它的名字不是WC和C，而是WC1，WC2，C1，C2，你的微分方程只需要W和C作为状态变量。

state <- c(WC = 0.5, C = 20)

这将给出您预期的结果！

顺便说一句:当你把所有的参数都放在一个向量"parms“中时，它会更清晰。

parms <- c(
   bulk = 0.5,
   FC = 0.4,
   Ks = 0.03,
   Sat = 0.8,
   Wres = 0.1,
   kbio = 0.01,
   kd = 10)

然后使用：

 out <- as.data.frame(ode(times=times, y=state, func=SPEC, parms=parms, method="rk4"))

祝你好运，约翰尼斯