Stan Math Library  2.20.0
reverse mode automatic differentiation
integrate_ode_rk45.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_ARR_FUNCTOR_INTEGRATE_ODE_RK45_HPP
2 #define STAN_MATH_PRIM_ARR_FUNCTOR_INTEGRATE_ODE_RK45_HPP
3 
13 #include <boost/version.hpp>
14 #if BOOST_VERSION == 106400
15 #include <boost/serialization/array_wrapper.hpp>
16 #endif
17 #include <boost/numeric/odeint.hpp>
18 #include <algorithm>
19 #include <ostream>
20 #include <functional>
21 #include <iterator>
22 #include <vector>
23 
24 namespace stan {
25 namespace math {
26 
73 template <typename F, typename T1, typename T2, typename T_t0, typename T_ts>
74 std::vector<std::vector<typename stan::return_type<T1, T2, T_t0, T_ts>::type>>
75 integrate_ode_rk45(const F& f, const std::vector<T1>& y0, const T_t0& t0,
76  const std::vector<T_ts>& ts, const std::vector<T2>& theta,
77  const std::vector<double>& x, const std::vector<int>& x_int,
78  std::ostream* msgs = nullptr,
79  double relative_tolerance = 1e-6,
80  double absolute_tolerance = 1e-6, int max_num_steps = 1E6) {
81  using boost::numeric::odeint::integrate_times;
82  using boost::numeric::odeint::make_dense_output;
83  using boost::numeric::odeint::max_step_checker;
84  using boost::numeric::odeint::runge_kutta_dopri5;
85 
86  const double t0_dbl = value_of(t0);
87  const std::vector<double> ts_dbl = value_of(ts);
88 
89  check_finite("integrate_ode_rk45", "initial state", y0);
90  check_finite("integrate_ode_rk45", "initial time", t0_dbl);
91  check_finite("integrate_ode_rk45", "times", ts_dbl);
92  check_finite("integrate_ode_rk45", "parameter vector", theta);
93  check_finite("integrate_ode_rk45", "continuous data", x);
94 
95  check_nonzero_size("integrate_ode_rk45", "initial state", y0);
96  check_nonzero_size("integrate_ode_rk45", "times", ts_dbl);
97  check_ordered("integrate_ode_rk45", "times", ts_dbl);
98  check_less("integrate_ode_rk45", "initial time", t0_dbl, ts_dbl[0]);
99 
100  if (relative_tolerance <= 0)
101  invalid_argument("integrate_ode_rk45", "relative_tolerance,",
102  relative_tolerance, "", ", must be greater than 0");
103  if (absolute_tolerance <= 0)
104  invalid_argument("integrate_ode_rk45", "absolute_tolerance,",
105  absolute_tolerance, "", ", must be greater than 0");
106  if (max_num_steps <= 0)
107  invalid_argument("integrate_ode_rk45", "max_num_steps,", max_num_steps, "",
108  ", must be greater than 0");
109 
110  // creates basic or coupled system by template specializations
111  coupled_ode_system<F, T1, T2> coupled_system(f, y0, theta, x, x_int, msgs);
112 
113  // first time in the vector must be time of initial state
114  std::vector<double> ts_vec(ts.size() + 1);
115  ts_vec[0] = t0_dbl;
116  std::copy(ts_dbl.begin(), ts_dbl.end(), ts_vec.begin() + 1);
117 
118  std::vector<std::vector<typename stan::return_type<T1, T2, T_t0, T_ts>::type>>
119  y;
120  coupled_ode_observer<F, T1, T2, T_t0, T_ts> observer(f, y0, theta, t0, ts, x,
121  x_int, msgs, y);
122  bool observer_initial_recorded = false;
123 
124  // avoid recording of the initial state which is included by the
125  // conventions of odeint in the output
126  auto filtered_observer
127  = [&](const std::vector<double>& coupled_state, double t) -> void {
128  if (!observer_initial_recorded) {
129  observer_initial_recorded = true;
130  return;
131  }
132  observer(coupled_state, t);
133  };
134 
135  // the coupled system creates the coupled initial state
136  std::vector<double> initial_coupled_state = coupled_system.initial_state();
137 
138  const double step_size = 0.1;
139  integrate_times(
140  make_dense_output(absolute_tolerance, relative_tolerance,
141  runge_kutta_dopri5<std::vector<double>, double,
142  std::vector<double>, double>()),
143  std::ref(coupled_system), initial_coupled_state, std::begin(ts_vec),
144  std::end(ts_vec), step_size, filtered_observer,
145  max_step_checker(max_num_steps));
146 
147  return y;
148 }
149 
150 } // namespace math
151 
152 } // namespace stan
153 
154 #endif
void check_finite(const char *function, const char *name, const T_y &y)
Check if y is finite.
void check_nonzero_size(const char *function, const char *name, const T_y &y)
Check if the specified matrix/vector is of non-zero size.
void check_ordered(const char *function, const char *name, const std::vector< T_y > &y)
Check if the specified vector is sorted into strictly increasing order.
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition: value_of.hpp:17
Observer for the coupled states.
std::vector< std::vector< typename stan::return_type< T1, T2, T_t0, T_ts >::type > > integrate_ode_rk45(const F &f, const std::vector< T1 > &y0, const T_t0 &t0, const std::vector< T_ts > &ts, const std::vector< T2 > &theta, const std::vector< double > &x, const std::vector< int > &x_int, std::ostream *msgs=nullptr, double relative_tolerance=1e-6, double absolute_tolerance=1e-6, int max_num_steps=1E6)
Return the solutions for the specified system of ordinary differential equations given the specified ...
void invalid_argument(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
Throw an invalid_argument exception with a consistently formatted message.
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:87
The coupled_ode_system represents the coupled ode system, which is the base ode and the sensitivities...
void check_less(const char *function, const char *name, const T_y &y, const T_high &high)
Check if y is strictly less than high.
Definition: check_less.hpp:63
const kernel_cl< in_buffer, out_buffer, int, int > copy("copy", {indexing_helpers, copy_kernel_code})
See the docs for copy() .

     [ Stan Home Page ] © 2011–2018, Stan Development Team.