r - コンパイル済みモデルの ODE システムの初期 (状態) 値の設定 (deSolve、Rcpp)

Question

コンパイルされた ODE を呼び出して R パッケージ'deSolve'を介して解決する際に、おそらくマイナーな問題に苦しんでおり、より専門的なユーザーからのアドバイスを求めています。

バックグラウンド

で解かなければならない ODE 系がいくつかあります'deSolve'。個別の C++ 関数 (モデルごとに 1 つ) で ODE を定義しました'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)
*/

出力は次のとおりです (ここでは最初の 2 行のみ)。

  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 が存在しなくても動作するはずですが、Windows では Rtools、Linux では GNU コンパイラ、macOS では Xcode が必要なコードをコンパイルするには、パッケージ'Rcpp'とパッケージをインストールする必要があります。'deSolve'

問題

私の理解では、(または) でneある必要があります。残念ながら、ソルバーは、ODE を除いて、C++ 関数で提供したものを考慮していないようです。ただし、R を別の値に変更すると、機能します。どういうわけか、C++ で定義した if 条件は無視されますが、R ではなく C++ で初期値を計算できる理由と方法がわかりません。0time = 1t < 2state