Ваше уравнение в B(t) почти разделимо, так как вы можете разделить B(t), из которого вы можете получить это

B(t) = C * exp{-p5 * t} * (p2 + B(t)) ^ {of_interest * p1 * p3}

Это неявное решение для B(t), которое мы будем решать точечно.

Вы можете решить для C, учитывая ваше начальное значение B. Я полагаю t = 0 изначально? В таком случае

C = B_0 / (p2 + B_0) ^ {of_interest * p1 * p3}

Это также дает несколько более красивое выражение для A(t):

dA(t) / dt = B_0 / (p2 + B_0) * p1 * p3 * (1 - of_interest) *
   exp{-p5 * t} * ((p2 + B(t) / (p2 + B_0)) ^ 
   {of_interest * p1 * p3 - 1} - p4 * A(t)

Это может быть решено путем интегрирования коэффициента (= exp{p4 * t}) путем численного интегрирования члена, включающего B(t). Мы указываем нижнюю границу интеграла как 0, чтобы нам никогда не приходилось вычислять B за пределами диапазона [0, t], что означает, что константа интегрирования просто A_0 и, таким образом:

A(t) = (A_0 + integral_0^t { f(tau; parameters) d tau}) * exp{-p4 * t}

Основная суть в том, что B(t) управляет всем в этой системе - подход будет следующий: определить поведение B(t), затем использовать это, чтобы выяснить, что происходит с A(t), а затем развернуть.

Во-первых, «внешние» параметры; нам также нужно nleqslv, чтобы получить B:

library(nleqslv)

t_min <- 0
t_max <- 10000
t_N <- 10

#we'll only solve the behavior of A & B over t_rng
t_rng <- seq(t_min, t_max, length.out = t_N)

#I'm calling of_interest ttheta
ttheta_min <- 0
ttheta_max <- 1
ttheta_N <- 5

tthetas <- seq(ttheta_min, ttheta_max, length.out = ttheta_N)

B_0 <- 1.4
A_0 <- 28

#No sense storing this as a vector when we'll only ever use it as a list
parameters <- list(p1 = 0.028, p2 = 0.3, p3 = 0.5, 
                   p4 = 0.0002, p5 = 0.001)

Отсюда основная схема:

1. Учитывая значения параметров (в частности, ttheta), найдите BB по t_rng с помощью решения нелинейного уравнения.
2. Учитывая BB и значения параметров, найдите AA больше t_rng путем численного интегрирования.
3. Учитывая AA и ваше выражение для dAdt, подключи и максимизируй.

производные ‹- sapply (tthetas, function (th) {#append current ttheta params‹ - c (parameters, ttheta = th)

#declare a function we'll use to solve for B (see above)
b_slv <- function(b, t) 
  with(params, b - B_0 * ((p2 + b)/(p2 + B_0)) ^ 
         (ttheta * p1 * p3) * exp(-p5 * t))

#solving point-wise (this is pretty fast)
#  **See below for a note**
BB <- sapply(t_rng, function(t) nleqslv(B_0, function(b) b_slv(b, t))$x)

#this is f(tau; params) that I mentioned above;
#  we have to do linear interpolation since the
#  numerical integrator isn't constrained to the grid.
#  **See below for note**
a_int <- function(t){
  #approximate t to the grid (t_rng)
  #  (assumes B is monotonic, which seems to be true)
  #  (also, if t ends up negative, just assign t_rng[1])
  t_n <- max(1L, which.max(t_rng - t >= 0) - 1L)
  idx <- t_n:(t_n+1)
  ts <- t_rng[idx]

  #distance-weighted average of the local B values
  B_app <- sum((-1) ^ (0:1) * (t - ts) / diff(ts) * BB[idx])
  #finally, f(tau; params)
  with(params, (1 - ttheta) * p1 * p3 * B_0 / (p2 + B_0) *
         ((p2 + B_app)/(p2 + B_0)) ^ (ttheta * p1 * p3 - 1) *
         exp((p4 - p5) * t))
}

#a_int only works on scalars; the numeric integrator
#  requires a version that works on vectors
a_int_v <- function(t) sapply(t, a_int)

AA <- exp(-params$p4 * t_rng) * 
  sapply(t_rng, function(tt)
    #I found the subdivisions constraint binding in some cases
    #  at the default value; no trouble at 1000.
    A_0 + integrate(a_int_v, 0, tt, subdivisions = 1000L)$value)

#using the explicit version of dAdt given as flux1 - flux2
max(with(params, (1 - ttheta) * p1 * p3 * BB / (p2 + BB) - p4 * AA))})

Finally, simply run `tthetas[which.max(derivs)]` to get the maximizer.

Примечание:

Этот код не оптимизирован для повышения эффективности. Есть несколько мест, где есть некоторые потенциальные возможности для ускорения:

• вероятно, быстрее рекурсивно запустить решатель уравнений, так как он будет быстрее сходиться с лучшими начальными предположениями - использование предыдущего значения вместо начального, безусловно, лучше
• Будет быстрее просто использовать суммы Римана для интеграции; компромисс заключается в точности, но все будет в порядке, если у вас достаточно плотная сетка. Одно из преимуществ Римана в том, что вам вообще не придется интерполировать, а численно это простая линейная алгебра. Я запускал это с помощью t_N == ttheta_N == 1000L, и он работал в течение нескольких минут.
• Вероятно, можно будет векторизовать a_int напрямую, а не просто sapply на нем, что сопутствует ускорению за счет более прямого обращения к BLAS.
• Куча другой мелочи. Выполните предварительное вычисление ttheta * p1 * p3, поскольку он так часто используется повторно и т. Д.

Однако я не стал включать что-либо из этого, потому что, честно говоря, вам, вероятно, лучше перенести это на более быстрый язык - Джулия - мой любимый питомец, но, конечно, R хорошо говорит с C ++, C, Fortran и т. Д. .