在编译模型 (deSolve, Rcpp) 中设置 ODE 系统的初始（状态）值

Question

在调用已编译的 ODE 来解决时，我正在努力解决一个可能很小的问题 通过 R 包'deSolve' 和我寻求更多专家用户的建议。

背景

我有几个 ODE 系统需要用 'deSolve' 解决。我已经在单独的 C++ 函数（每个模型一个）中定义了 ODE，我通过 R 与 'Rcpp' 一起调用。如果函数从另一个模型中获取输入，系统的初始值会发生变化（因此基本上是级联）。

这很好用，但是，对于一个模型，我必须为t < 2 设置初始参数。我尝试在 C++ 函数中执行此操作，但它似乎不起作用。

运行代码示例

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export("set_ODE")]]
SEXP set_ODE(double t, NumericVector state, NumericVector parameters) {
  
  List dn(3);
  
  double tau2 = parameters["tau2"]; 
  double Ae2_4 = parameters["Ae2_4"]; 
  double d2 = parameters["d2"]; 
  double N2 = parameters["N2"];
  
  double n2 = state["n2"];
  double m4 = state["m4"];
  double ne = state["ne"];
  
  // change starting conditions for t < 2
  if(t < 2) {
    n2 = (n2 * m4) / N2;
    m4 = n2;
    ne = 0;
    
  }
  
  dn[0] = n2*d2 - ne*Ae2_4 - ne/tau2;
  dn[1] = ne/tau2 - n2*d2;
  dn[2] = -ne*Ae2_4;    
  
  return(Rcpp::List::create(dn));
}


/*** R
state <-  c(ne = 10, n2 = 0, m4 = 0)
parameters <- c(N2 = 5e17, tau2 = 1e-8, Ae2_4 = 5e3, d2 = 0)

results <- deSolve::lsoda(
  y = state,
  times = 1:10,
  func = set_ODE,
  parms = parameters
)

print(results)
*/

输出读取（这里只有前两行）：

  time            ne           n2            m4
1     1  1.000000e+01 0.000000e+00  0.000000e+00
2     2  1.000000e+01 2.169236e-07 -1.084618e-11

以防万一：如何运行此代码示例？

我的示例使用RStudio进行了测试：

• 将代码复制到以 *.cpp 结尾的文件中
• 点击“来源”按钮（或<shift> + <cmd> + <s>）

它应该在没有 RStudio 存在的情况下也能工作，但必须安装包 'Rcpp' 和 'deSolve' 并编译它需要在 Windows 上使用 Rtools、在 Linux 上使用 GNU 编译器和在 macOS 上使用 Xcode 的代码。

问题

据我了解，ne 对于time = 1（或t < 2）应该是0。不幸的是，求解器似乎没有考虑我在 C++ 函数中提供的内容，除了 ODE。但是，如果我将 R 中的 state 更改为另一个值，它就可以工作。不知何故，我在 C++ 中定义的 if 条件被忽略了，但我不明白为什么以及如何在 C++ 而不是 R 中计算初始值。

Answer 1

我能够重现您的代码。在我看来，这确实很优雅，即使它没有利用求解器的全部功能。原因是，Rcpp 通过一个普通的 R 函数为编译模型创建了一个接口。因此，在每个时间步骤中，从 slover（例如 lsoda）到 R 的回拨是必要的。这种反向调用不适用于“普通”C/Fortran 接口。这里求解器和模型之间的通信发生在机器代码级别。

通过这些信息，我可以看到我们不需要在 C/C++ 级别上期待初始化问题，但它看起来像是一个典型案例。由于模型函数只是模型的导数（并且仅此而已）。积分由求解器“从外部”完成。它始终使用从之前的时间步（粗略地说）派生的实际积分状态调用模型。因此，无法在模型函数中强制状态变量为固定值。

但是，有几种方法可以解决这个问题：

• lsoda 调用链接
• 使用事件

下面展示了一个链式的方法，但是我还不确定第一个时间段的参数的初始化，所以可能只是解决方案的一部分。

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export("set_ODE")]]
SEXP set_ODE(double t, NumericVector state, NumericVector parameters) {

  List dn(3);

  double tau2 = parameters["tau2"];
  double Ae2_4 = parameters["Ae2_4"];
  double d2 = parameters["d2"];
  double N2 = parameters["N2"];

  double n2 = state["n2"];
  double m4 = state["m4"];
  double ne = state["ne"];

  dn[0] = n2*d2 - ne*Ae2_4 - ne/tau2;
  dn[1] = ne/tau2 - n2*d2;
  dn[2] = -ne*Ae2_4;

  return(Rcpp::List::create(dn));
}


/*** R
state <-  c(ne = 10, n2 = 0, m4 = 0)
parameters <- c(N2 = 5e17, tau2 = 1e-8, Ae2_4 = 5e3, d2 = 0)

## the following is not yet clear to me !!!
## especially as it is essentially zero
y1 <- c(ne = 0,
       n2 = unname(state["n2"] * state["m4"]/parameters["N2"]),
       m4 = unname(state["n2"]))


results1 <- deSolve::lsoda(
  y = y,
  times = 1:2,
  func = set_ODE,
  parms = parameters
)

## last time step, except "time" column
y2 <- results1[nrow(results1), -1]

results2 <- deSolve::lsoda(
  y = y2,
  times = 2:10,
  func = set_ODE,
  parms = parameters
)

## omit 1st time step in results2
results <- rbind(results1, results2[-1, ])

print(results)
*/

该代码还有另一个潜在问题，因为参数跨越从 1e-8 到 1e17 的几个数量级。这可能会导致数值问题，因为包括 R 在内的大多数软件的相对精度仅涵盖 16 个数量级。这可能是原因，为什么结果都是零？在这里重新调整模型可能会有所帮助。