求解 R 中的微分方程组

Question

我在 R 中有一个简单的通量模型。它归结为两个微分方程，对模型中的两个状态变量建模，我们将它们称为 A 和 B。它们被计算为四个分量通量 flux1-flux4、5 个参数 p1-p5 和第 6 个参数 of_interest 的简单差分方程，其值可以介于 0-1 之间。

parameters<- c(p1=0.028, p2=0.3, p3=0.5, p4=0.0002, p5=0.001, of_interest=0.1) 
state     <- c(A=28, B=1.4)

model<-function(t,state,parameters){
  with(as.list(c(state,parameters)),{
  #fluxes
  flux1  = (1-of_interest) * p1*(B / (p2 + B))*p3
  flux2  = p4* A          #microbial death
  flux3  = of_interest * p1*(B / (p2 + B))*p3 
  flux4  = p5* B      

  #differential equations of component fluxes
  dAdt<- flux1 - flux2
  dBdt<- flux3 - flux4
  list(c(dAdt,dBdt))
  })

我想写一个函数来求dAdt关于of_interest的导数，将导出的方程设为0，然后重新排列并求解[=17=的值].这将是最大化函数 dAdt.

的参数 of_interest 的值

到目前为止，我已经能够在稳定状态下求解模型，跨 of_interest 的可能值来证明应该有一个最大值。

require(rootSolve)
range<- seq(0,1,by=0.01)
for(i in range){
of_interest=i
parameters<- c(p1=0.028, p2=0.3, p3=0.5, p4=0.0002, p5=0.001, of_interest=of_interest) 
state     <- c(A=28, B=1.4)
ST<- stode(y=y,func=model,parms=parameters,pos=T)
out<- c(out,ST$y[1])

然后绘图：

plot(out~range, pch=16,col='purple')
lines(smooth.spline(out~range,spar=0.35), lwd=3,lty=1)

我如何分析求解 of_interest 的值使 R 中的 dAdt 最大化？如果解析解是不可能的，我怎么知道，我怎样才能用数值方法解决这个问题？

更新：我认为这个问题可以通过 R 中的 deSolve 包解决，链接 here，但是我在使用我的特定示例实现它时遇到了问题。

Answer 1

你在 B(t) 中的等式是 just-about 可分的，因为你可以除掉 B(t)，从中你可以得到

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

这是 B(t) 的隐式解决方案，我们将解决 point-wise。

您可以根据 B 的初始值求解 C。我想 t = 0 最初？在这种情况下

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

这也为 A(t) 提供了一个有点 nicer-looking 的表达式：

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)

这可以通过对涉及 B(t) 的项进行数值积分，通过积分因子 (= exp{p4 * t}) 来解决。我们将积分的下限指定为 0，这样我们就不必在运行ge [0, t] 之外计算 B，这意味着积分常数只是 A_0，因此：

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

基本要点是 B(t) 正在驱动这个系统中的一切——方法是：解决 B(t) 的行为，然后用它来弄清楚 [= 发生了什么22=]，然后最大化。

首先，"outer"参数；我们还需要 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)

从这里开始，基本轮廓是：

给定参数值（特别是 ttheta），通过 non-linear 方程求解

BB

t_rng

给定 BB 和参数值，通过数值积分 t_rng 求解 AA
给定 AA 和你的 dAdt 表达式，插入并最大化。

派生 <- 应用（tthetas，函数（th）{ #append 当前 ttheta 参数 <- c（参数，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.

注：

此代码未针对效率进行优化。有几个地方有一些潜力speed-ups:

递归地运行方程求解器可能更快，因为它会更快地收敛并提供更好的初始猜测——使用先前的值而不是初始值肯定更好
单纯使用黎曼和积分会更快；权衡是在准确性上，但如果你有一个足够密集的网格应该没问题。黎曼的一个优点是你根本不需要插值，而且在数值上它们是简单的线性代数。我运行使用 t_N == ttheta_N == 1000L 并在几分钟内完成运行。
可能可以直接矢量化 a_int 而不是仅 sapply 对其进行矢量化，同时 speed-up 通过更直接地诉诸 BLAS。
大量其他小东西。 Pre-compute ttheta * p1 * p3 因为re-used 这么多，等等

不过，我没有费心包括任何这些东西，因为老实说，将它移植到一种更快的语言可能会更好——Julia 是我自己最喜欢的，但当然 R 与 C++ 交流得很好， C、要塞运行等

求解 R 中的微分方程组

Solving a system of differential equations in R

r

optimization

differential-equations

calculus

maximization

注：