Как получить результаты с фиксированным шагом?

Question

У меня есть система дифференциальных уравнений, она немного сложна, но самая близкая аналогия, к которой я могу прийти, - это каталитическая химическая реакция, при которой катализатор разлагается со временем, например

      A + B -> C + B  (rate k1)
B -> D          (rate k2)

Итак, у нас есть а также dBbt(a,B) = -k2*B

Я могу смоделировать это с помощью и фиксированный набор временных точек ... но я хочу запустить это, пока не скажу B <1e-8

Я думаю, что это можно сделать, используя корневые функции в общем пакет, но это также требует только параметр, который заставляет меня заранее выбирать время, в течение которого я хочу «концентрации A и B».

Примечание: мои настоящие уравнения более сложны и на самом деле их довольно сложно вычислить - они не имеют ничего общего с химией, кроме значений, все являются положительными числами, и как только одно из них становится настолько близко к нулю, что не имеет никаких шансов, состояние системы перестает меняться. .

У меня есть скрученная вручную реализация RK4, которая делает то, что я хочу, но код, который я пишу, - это код, который мне нужно протестировать, я знаю, что пакет хорошо протестирован, поэтому я просто ищу способ получить выходные данные с фиксированными размерами шагов до тех пор, пока "корневая" функция возвращает истину

Обновлять

Вот пример решения моей проблемы, когда я знаю, когда мне нужен ответ:

      kernel <- function(t, y, params, input) {
  with(as.list(c(params,y)), {
    dA <- -k1*A*B
    dB <- -k2*B
    list(c(dA,dB))
  })
}
params <- c(k1=0.03, k2=0.005)
y0 <- c(A=1.3, B=0.5)

# here I need to pick the times!
times <- seq(0,500, by=1)
out <- rk4(y0, times, kernel, params)

Я хочу что-то вроде «Замените последние две строки кода на это».

      out <- ___name_of_function___(
  y0, 
  initial_time=0, 
  delta_time=1, 
  kernel, 
  params, 
  rootfun=function(t,y,p){y.B<1e-8}
)

Где надеюсь какая-то функция в пакет или связанный вспомогательный пакет

Answer 1

Пример 2 в https://www.rdocumentation.org/packages/deSolve/versions/1.28/topics/lsodar показывает, что общий подход намекает, что для двух вызовов решателя мы можем получить второй вызов для оценки всех требуемых точек, за счет отсутствия контроля над количеством оценок при первом обращении:

      kernel <- function(t, y, params, input) {
  with(as.list(c(params,y)), {
    dA <- -k1*A*B
    dB <- -k2*B
    list(c(dA,dB))
  })
}
params <- c(k1=0.03, k2=0.005)
y0 <- c(A=1.3, B=0.5)
rootfun <- function(t, y, params, input) {
  with(as.list(c(params,y)), {
    ifelse(B<1e-8,0,B)
  })
}

# First find where the root is
out <- ode(y0, c(0,1e8), kernel, params, root=rootfun)
# Now the second rwo of `out` will be the location of the root
# So just create the times up to that point
times <- seq(0,round(out[2,'time']), by=1)
out <- ode(y0, times, kernel, params, root = rootfun)

ПРИМЕЧАНИЕ: количество оценок функции во время поиска корня ограничено maxsteps, которое по умолчанию равно 500.