使用Rcpp求解具有时变参数的ODE

Question

我的目标是使用Rcpp求解微分方程组。基本上我想建立一个系统，如下面的代码所示（修改代码示例：How to use C++ ODE solver with Rcpp in R?）。

目前，下面的代码在时间间隔0到10中整合了一组颂歌。整个时间段params [0]是-100，而parms [1] = 10.但是，我的目标是设置一个系统parms [0]和parms [1]仅在时间间隔的子集上保持不变。例如。对于时间间隔0-5 parms [0]应设置为1，剩余时间parms [0]应为10。

实际上，我几乎没有c ++ / rcpp的经验。因此，我不知道如何建立这样一个系统。你能不能给我一个暗示我应该如何构建颂歌系统的呢？非常感谢您提出如何解决此问题的任何建议。

我将下面的代码保存在cpp文件中并使用R中的sourceCpp调用它：

#include <Rcpp.h>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>

// [[Rcpp::depends(BH)]]
using namespace Rcpp;
using namespace std;
using namespace boost::numeric::odeint;

typedef boost::array< double ,3 > state_type;
typedef boost::array< double ,2 > parms_type;

double time = 10;
parms_type parms = {-100, 10};

void rhs( const state_type &x , state_type &dxdt , const double t) {
dxdt[0] = parms[0]/(2.0*t*t) + x[0]/(2.0*t);
dxdt[1] = parms[1]/(2.0*t*t) + x[1]/(2.0*t);
dxdt[2] = parms[1]/(2.0*t*t) + x[1]/(2.0*t);
}

void write_cout( const state_type &x , const double t ) {
// use Rcpp's stream
Rcpp::Rcout << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2] <<  endl;
}

typedef runge_kutta_dopri5< state_type > stepper_type;

// [[Rcpp::export]]
bool boostExample() {
state_type x = { 1.0 , 1.0, 1.0 }; // initial conditions
integrate_adaptive(make_controlled( 1E-12 , 1E-12 , stepper_type () ) ,
                 rhs , x , 1.0 , time, 0.1 , write_cout );
return true;
}

Answer 1

你的代码不能为我编译：

boost-ode.cpp:11:8: error: ‘double time’ redeclared as different kind of symbol
 double time = 10.0;
        ^~~~
In file included from /usr/include/pthread.h:24:0,
                 from /usr/include/x86_64-linux-gnu/c++/6/bits/gthr-default.h:35,
                 from /usr/include/x86_64-linux-gnu/c++/6/bits/gthr.h:148,
                 from /usr/include/c++/6/ext/atomicity.h:35,
                 from /usr/include/c++/6/bits/basic_string.h:39,
                 from /usr/include/c++/6/string:52,
                 from /usr/include/c++/6/stdexcept:39,
                 from /usr/include/c++/6/array:39,
                 from /usr/include/c++/6/tuple:39,
                 from /usr/include/c++/6/unordered_map:41,
                 from /usr/local/lib/R/site-library/Rcpp/include/Rcpp/platform/compiler.h:153,
                 from /usr/local/lib/R/site-library/Rcpp/include/Rcpp/r/headers.h:48,
                 from /usr/local/lib/R/site-library/Rcpp/include/RcppCommon.h:29,
                 from /usr/local/lib/R/site-library/Rcpp/include/Rcpp.h:27,
                 from boost-ode.cpp:1:
/usr/include/time.h:192:15: note: previous declaration ‘time_t time(time_t*)’
 extern time_t time (time_t *__timer) __THROW;
               ^~~~

我只是删除了全局变量time并在其中使用了明确的10.0。我还删除了Rcpp和std的命名空间用法。无论如何都没有使用前者，后者只在一个地方使用。通常我会尽量避免导入这么大的命名空间，特别是两个同时。

无论如何，一个简单的解决方案是引入两个参数向量，并根据时间在rhs中选择合适的一个：

#include <Rcpp.h>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>

// [[Rcpp::depends(BH)]]
using namespace boost::numeric::odeint;

typedef boost::array< double ,3 > state_type;
typedef boost::array< double ,2 > parms_type;

parms_type parms_begin = {1, 10};
parms_type parms_end = {10, 10};

void rhs( const state_type &x , state_type &dxdt , const double t) {
  if (t < 5.0) {
    dxdt[0] = parms_begin[0]/(2.0*t*t) + x[0]/(2.0*t);
    dxdt[1] = parms_begin[1]/(2.0*t*t) + x[1]/(2.0*t);
    dxdt[2] = parms_begin[1]/(2.0*t*t) + x[1]/(2.0*t);
  } else {
    dxdt[0] = parms_end[0]/(2.0*t*t) + x[0]/(2.0*t);
    dxdt[1] = parms_end[1]/(2.0*t*t) + x[1]/(2.0*t);
    dxdt[2] = parms_end[1]/(2.0*t*t) + x[1]/(2.0*t);
  }
}

void write_cout( const state_type &x , const double t ) {
  // use Rcpp's stream
  Rcpp::Rcout << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2] <<  std::endl;
}

typedef runge_kutta_dopri5< state_type > stepper_type;

// [[Rcpp::export]]
bool boostExample() {
  state_type x = { 1.0 , 1.0, 1.0 }; // initial conditions
  integrate_adaptive(make_controlled( 1E-12 , 1E-12 , stepper_type () ) ,
                     rhs , x , 1.0 , 10.0, 0.1 , write_cout );
  return true;
}
/*** R
boostExample()
*/