在编译模型（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结尾的文件中
单击 'Source' 按钮（或 <shift> + <cmd> + <s>）

它也应该在没有 RStudio 的情况下工作，但是必须安装包 'Rcpp' 和 'deSolve' 并且编译它需要的代码 Windows 上的 Rtools，[= 上的 GNU 编译器58=] 和 Xcode 在 macOS 上。

问题

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

Answer 1

我能够重现您的代码。在我看来，这确实很优雅，即使它没有充分利用求解器的全部功能。原因是，Rcpp 通过普通的 R 函数创建了一个已编译模型的接口。因此，在每个时间步长中，从 slovers（例如 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 个数量级。难道这就是原因，为什么结果都是零？这里可能有助于重新缩放模型。

在编译模型（deSolve、Rcpp）中为 ODE 系统设置初始（状态）值

Setting initial (state) values for ODE system in compiled model (deSolve, Rcpp)

r

desolve

背景

运行代码示例

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

问题

在编译模型（deSolve、Rcpp）中为 ODE 系统设置初始（状态）值

Setting initial (state) values for ODE system in compiled model (deSolve, Rcpp)

r

desolve

背景

运行 代码示例

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

问题

运行代码示例

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