Я пытаюсь решить дифференциальное уравнение первого порядка, используя функцию ode из пакета deSolve. Проблема заключается в следующем: лекарство вводят с постоянной скоростью инфузии в несколько раз (время инфузии) и устраняют в скорости первого порядка. Таким образом, процесс может быть описан следующим образом:
if(t %in% Infusion_times){Infusion <- Infusion_rate} else{Infusion <- 0}
dC <- -Ke*C + Infusion
, где t - время, Infusion_times - вектор, содержащий время, в которое препарат вводится, C - количество препарата. , Ke - это ее постоянная элиминации, а Infusion - это переменная, принимающая значение скорости инфузии, когда есть инфузия, и значение 0 в противном случае. Итак, допустим, мы хотим назначить 3 дозы, начиная с 0, 24 и 40, при этом каждая инфузия длится два часа, и, скажем, мы хотим, чтобы deSolve вычислял ответ каждые 0,02 единицы времени. Мы хотим, чтобы deSolve решил дифференциальное уравнение для времен от 0 до 48, например, с шагом 0,02. Таким образом, наш вектор, определяющий время для функции ode, будет иметь вид:
times <- seq(from = 0, to = 48, by = 0.02)
Время инфузии определяется следующим образом:
Infusion_times <- c(seq(from = 0, to = 2, by = 0.02), seq(from = 24, to = 26, by = 0.02),
seq(from = 40, to = 42, by = 0.02))
Сначала я был обеспокоен тем, что проблема может быть в сравнении с плавающей точкой. Чтобы предотвратить это, я округлил оба вектора до двух десятичных разрядов:
times <- round(times, 2)
Infusion_times <- round(times, 2)
Так что теперь, надеюсь, все Infusion_times включены в вектор times:
(sum(Infusion_times %in% times)/length(Infusion_times))*100
[1] 100
As Вы можете видеть, что все значения в Infusion_times (100%) содержатся в векторе times, и поэтому переменная Infusion должна принимать значение Infusion_rate в указанное время. Однако, когда мы решаем уравнение, оно не работает. Давайте докажем это, но сначала давайте определим другие значения, необходимые для функции ode:
parameters <- c(Ke = 0.5)
amounts <- c(C = 0) #Initial value for drug is 0
Inf_rate <- 5
А теперь давайте напишем функцию, определяющую скорость изменений, так же, как требуется:
OneComp <- function(t, amounts, parameters){
with(as.list(c(amounts, parameters)),{
if(t %in% Infusion_times){Infuse =Inf_rate} else{Infuse = 0}
dC <- -Ke*C + Infuse
list(c(dC))})
}
Для тех, кто не знаком с пакетом deSolve, аргумент t функции будет указывать время, за которое должно быть интегрировано уравнение, amounts будет указывать начальное значение C и parameters даст значение параметров (в данном случае просто Ke). А теперь давайте решим уравнение:
out <- ode(func = OneComp, y = amounts, parms = parameters, times = times)
Давайте построим результаты:
library(ggplot2)
ggplot(as.data.frame(out)) + geom_line(aes(x = times, y = C))
Это точно такой же результат, который мы получили бы, если бы Infusion всегда было равно 0. Тем не менее, обратите внимание, что если бы мы имитировали c только одну дозу, и мы должны были попробовать подобный подход, он бы работал:
OneComp <- function(t, amounts, parameters){
with(as.list(c(amounts, parameters)),{
if(t < 2){Infuse =Inf_rate} else{Infuse = 0}
dC <- -Ke*C + Infuse
list(c(dC))})
}
out <- ode(func = OneComp, y = amounts, parms = parameters, times = times)
ggplot(as.data.frame(out)) + geom_line(aes(x = times, y = C))
Здесь мы сделали переменную Infuse take значение Inf_rate только тогда, когда время меньше 2 часов, и оно работает! Поэтому я совершенно озадачен этим поведением. Это не проблема изменения значения переменной, это не проблема сравнения чисел с плавающей запятой ... Есть идеи, что это может быть? Спасибо