Решение ODE с изменяющимся во времени параметром с использованием Rcpp

Question

Моя цель - решить систему дифференциальных уравнений с использованием Rcpp. По сути, я хочу настроить систему, как показано в приведенном ниже коде (модификация примера кода приведена здесь: Как использовать C ++ ODE Solver с Rcpp в R? ).

В настоящий момент приведенный ниже код объединяет набор од в интервале времени от 0 до 10. Для всего времени параметры [0] составляют -100, а параметры [1] = 10. Однако моя цель - настроить система, в которой parms [0] и parms [1] постоянны только в подмножестве временного интервала. Например. для промежутка времени 0-5 parms [0] должно быть установлено значение 1, а для оставшегося времени parms [0] должно быть 10.

На самом деле, у меня почти нет опыта работы с c ++ / rcpp. Таким образом, я понятия не имею, как настроить такую ​​систему. Не могли бы вы дать мне подсказку, как я должен построить систему Оде. Заранее большое спасибо за любой совет, как решить эту проблему.

Я сохраняю код ниже в файле cpp и вызываю его с sourceCpp в R:

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