R中的deSolve包，变量不更新答案

【问题标题】：deSolve package in R, variable not updatingR中的deSolve包，变量不更新
【发布时间】：2016-04-28 16:35:00
【问题描述】：

我正在使用 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），求解器仅使用WC（0.5）的初始值，因此每次迭代计算相同的值（perc1和perc2始终等于0.007，这对于WC=0.5是正确的）。

奇怪的是，当我使用单层时，情况并非如此。此问题接缝类似于此处报告的内容： Population values not updating in deSolve in R 但是，我以“状态=值”的方式定义初始值，我仍然面临更新问题。有什么想法可以解决这个问题以及为什么我会面临这个问题？

【问题讨论】：

标签： r ode

【解决方案1】：

尝试将perc 的分配放在ifelse 之外。

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 中做作业会做什么。

【讨论】：

感谢您的推荐，我将分配移到 ifelse 之外。结果与我上面发布的相同：如果我使用单层，则所有内容都经过适当计算，2 层“perc”是仅使用初始 WC 值 0.5 计算

【解决方案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"))

最好的祝愿，约翰内斯

【讨论】：