您的方程式 inB(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)
从这里开始,基本轮廓是:
- 给定参数值(特别是
ttheta
),BB
通过t_rng
非线性方程求解求解
- 给定和参数值,通过数值积分
BB
求解AA
t_rng
- 给定
AA
和您对 dAdt 的表达,即插即用和最大化。
derivs <- 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
,因为它被重复使用了这么多,等等。
不过,我没有费心包括任何这些东西,因为老实说,你最好将它移植到更快的语言——Julia 是我自己的宠物,但当然 R 可以很好地与 C++、C、Fortran 等配合使用.