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

Question

我的目标是使用 Rcpp 求解微分方程组。基本上我想设置一个系统，如下面的代码所示（对此处找到的代码示例的修改：）。

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

其实我对c++/rcpp几乎没有经验。因此，我不知道如何建立这样一个系统。你能给我一个提示，我应该如何构建 ode 系统。非常感谢您提供有关如何解决此问题的任何建议。

我将下面的代码保存在一个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 中引入两个参数向量和 select，根据时间选择合适的向量：

#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()
*/

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

Solving ODE with time varying parameter using Rcpp

rcpp