пакет deSolve в R, переменная, не обновляющаяся

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 

Результаты для концентрации в двух слоях почвы (С1 и С2), являются правильными и идентичны результатам одного слоя. Вычисленное содержание воды, однако, неверно (результат для первого слоя должен быть идентичен симуляции с использованием одного слоя). Он швы, чтобы вычислить перколяцию (perc), решатель использует только начальное значение WC (0.5), поэтому одно и то же значение вычисляется на каждую итерацию (perc1 и perc2 всегда равны 0,007, что верно для WC = 0,5).

Странно, что это был не тот случай, когда я использовал один слой. Эта проблема швы похожа на то, что было сообщено здесь: Population values not updating in deSolve in R Однако я определяю начальные значения в форме «состояние = значение», и я все еще сталкиваюсь с проблемой обновления. Любые идеи, как я мог бы решить это и почему я столкнулся с этим?

Answer 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.

Answer 2

Ваша проблема в том, что вы неправильно определили свои переменные состояния. Выполнить свой код, а затем

state 

Вы видите, что имя не туалет и C, но WC1 WC2 C1 и C2 Ваше дифференциальное уравнение нужно только W и C в качестве переменных состояния.

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

Это дает ожидаемый результат!

BTW: Это понятнее, когда вы сложите все ваши параметры в векторе "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")) 

С наилучшими пожеланиями, Johannes