Эти проблемы могут быть сложными для диагностики, но я начал с рефакторинга вашей функции - то есть я упростил кишки и убедился, что он дал достаточно близкие ответы (используя all.equal()), чтобыоригинал.
- ваши первые 7 значений градиента используют идентичные формы и могут быть объединены в один векторизованный вызов.Это более эффективно, легче читать и легче отлаживать.
- , чтобы было легче увидеть, где возникает проблема, я далее разбил выражение на ряд более простых выражений (некоторые из которых повторяются висходный вызов: см. № 1 о преимуществах уменьшения дублирования ...)
- вычеркнул ненужное
tmp3^2/tmp3 (== tmp3) в уравнении - Я поставил
if(any(is.na(grad))) test и browser() вызывают функцию, чтобы мы могли остановиться при появлении первого значения NA/NaN и посмотреть, что происходит ...
func2 <- function(t,state,parameters, debug=TRUE){
n <- length(state)
v <- 1:(n-1)
grad <- rep(NA,n)
tmp1 <- (4*3.14*rho[v]*N0[v])
tmp2 <- 3*state[v]/tmp1
tmp3 <- tmp2^(1/3)
grad[v] <- with(as.list(c(state,parameters)),{
-D*(N0[v]*4*3.14*tmp3)*(Cs-S/V)
})
grad[n] <- -sum(grad[v])
if (debug && any(is.na(grad))) browser()
return(list(grad))
}
## test near-equality
all.equal(func1(0,state, parameters),func2(0,state, parameters)) ## TRUE
Теперь попробуйте запустить
out <- ode(y = state, times = times,func = func2, parms = parameters)
Это переводит нас в среду интерактивного браузера.
Первое промежуточное выражение выглядит нормально (большое, но конечное):
Browse[2]> tmp1
[1] 8724442 28341529 121926846 347177124 640918307 1295801866
[7] 11127053948
Второе выражение (3*state[v]/tmp1) выглядитХорошо, но обратите внимание последнее значение отрицательное - это, вероятно, потому что последняя (седьмая) переменная состояния стала слегка отрицательной.
Browse[2]> tmp2
X1 X2 X3 X4 X5
1.289771e-11 6.262837e-12 1.333549e-12 3.037421e-13 2.588684e-14
X6 X7
3.751315e-15 -4.992697e-18
Теперь, когда мы пытаемся получить корень кубадела идут плохо: если ценность неоштрафованный как тип complex, дробные степени отрицательных чисел равны NaN в R
Browse[2]> tmp3
X1 X2 X3 X4 X5 X6
2.345151e-04 1.843276e-04 1.100702e-04 6.722049e-05 2.958192e-05 1.553798e-05
X7
NaN
Это быстро распространится и испортит все состояние.
На этом этапе мымог бы попытаться отследить немного дальше и попытаться понять неточность с плавающей точкой, которая в первую очередь привела к отрицательному значению;однако может быть, а может и не быть легко (или даже возможно) переписать выражения таким образом, чтобы они были достаточно стабильными. В этом вопросе и в этом вопросе обсуждаются способы ограничения решений ОДУ неотрицательными значениями ... самое простое (если это имеет смысл для вашей проблемы) - это ввести pmax(tmp3,0) или pmax(tmp3,very_small_positive_number) вызов для предотвращения отрицательных значений ...
Незначительные комментарии:
- являются ли коэффициенты 3,14 в ваших вопросах предполагаемыми пи?R имеет встроенное
pi значение ... - для данного набора параметров, похоже, что
tmp1 является постоянным во времени.Возможно, вы захотите предварительно вычислить его вне функции градиента для эффективности ...
Чтобы увидеть, что происходит, я добавил na.rm=TRUE к сумме, как вы и предлагали.Я переключился на method="euler";это менее эффективно, но проще, поскольку делает очень мало промежуточных вычислений между вычислениями градиента.
out <- ode(y = state, times = times,func = func2, parms = parameters,
debug=FALSE,method="euler")
out <- out[rowSums(is.na(out))<9,]
png("SO_ode.png")
par(las=1)
matplot(out[,-1],type="l",lty=1,log="xy",col=1:8,
xlab="time",ylab="")
dev.off()
Это показывает, что один компонент за другимпадение до очень маленьких значений (и становится NaN после того, как мы пытаемся получить корень куба отрицательного значения ...) В зависимости от того, что вы делали, может быть безопасным для установки градиентов, значения которых былиNaN до нуля ...