2017-05-16 10 views
0

こんにちは私はこのODEを動作させるために戦ってきましたが、私はこのエラーを克服し続けます:eval(expr、envir、enclos)のエラー:オブジェクト 'j'が見つかりませんは微分方程式にif文を置くことができません

私のコードは以下の通りです、ODE

parameters <- c(
       a = 0.032, 
       b = (9/140), 
       c = (5/1400), 
       d = (95/700), 
       k = 1/140, 
       i = 0.25, 
       # r = 0.2, 
       n = 6000000, 
       x = 0.3 , 
       t = 1/180,  # important in looking at the shape 
       u = 1/180,  # important in looking at the shape 
       v = 1/360,  # important in looking at the shape 
       p = 10, 
       s = 10000, 
       g = 100 

       # e = .4, 
       #h = 1000 
) 

state <- c(
      S = 5989900, 
      E = 100, 
      I = 100, 
      Q = 100, 
      D = 100, 
      B = 100, 
      C = 100, 
      Y = 100, 
      H = 1000, 
      R = 1000, 
      J = 1000, 
      h = 100, 
      e = 0.1, 
      r = 0.1 
) 


# set up the equations 

equation <- (function(t, state, parameters) 
    with(as.list(c(state, parameters)), { 
    # rate of change 

    dS <- (-(a * S * I)/n) - (((1/r) * S * D)/n) 
    dE <- (a * S * I)/n + (((1/r) * S * D)/n) - i * E 
    if (h > Q) 
     j = 1 
    else if (h < Q) 
     j = 0 
    dI <- i * (j) * E - (e) * I - c * I - d * I 
    dQ <- (j) * (e) * I - b * Q - k * Q 
    dD <- d * I - r * D 
    dB <- b * Q + r * D 
    dC <- c * I + k * Q 

    dY <- p * (b * Q + r * D) 
    dR <- (1 - x) * (p * (b * Q + r * D)) - t * (R) 
    de <- t * (s/R) 
    dJ <- (x) * (p * (b * Q + r * D)) - v * (J) 
    dr <- v * (s/J) 
    dH <- (x) * (p * (b * Q + r * D)) - u * (H) 
    dh <- u * (H/g) 


    # return the rate of change 
    list(c(dS, dE, dI, dQ, dD, dB, dC, dY, dR, de, dJ, dr, dH, dh)) 
    })) 
# 

# solve the equations for certain starting parameters 


library(deSolve) 
times <- seq(0, 200, by = 1) 

out <- 
    ode(y = state, 
    times = times, 
    func = equation, 
    parms = parameters 
) 
# , method = "vode" 
head(out) 
tail(out) 

# graph the results 

par(oma = c(0, 0, 3, 0)) 
plot(out, xlab = "Time", ylab = "People") 
#plot(out[, "X"], out[, "Z"], pch = ".") 
mtext(outer = TRUE, side = 3, "Ebola Model",cex = 1.5 
) 

任意のヘルプは素晴らしいことだで私のif文の問題のようです!

+0

私は 'J' の宣言が表示されていないが、私は 'J' 参照してください。そして、宣言されていない場合は、文の中で宣言することはできません。 – Adamm

答えて

2

h==Q変数jが作成されない場合。 与えられた例では、hは、Qと等しくなります。 elseステートメントを追加するか、基本値をにifステートメントの前に割り当てる必要があります。このよう

:それは誤りでだのと同じように、

j = 0 

if (h > Q){ 
    j = 1 
} 
else if (h < Q) { 
    j = 0 
} 

または

if (h > Q){ 
    j = 1 
}else if (h < Q) { 
    j = 0 
}else{ 
j = 0 
} 
+0

を使うか、 '> ='または '<='(厳密ではない)の不等式 –