1 #ifndef STAN_MATH_PRIM_ARR_FUNCTOR_INTEGRATE_ODE_RK45_HPP 2 #define STAN_MATH_PRIM_ARR_FUNCTOR_INTEGRATE_ODE_RK45_HPP 13 #include <boost/version.hpp> 14 #if BOOST_VERSION == 106400 15 #include <boost/serialization/array_wrapper.hpp> 17 #include <boost/numeric/odeint.hpp> 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>>
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 = 1
e-6,
80 double absolute_tolerance = 1
e-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;
87 const std::vector<double> ts_dbl =
value_of(ts);
90 check_finite(
"integrate_ode_rk45",
"initial time", t0_dbl);
92 check_finite(
"integrate_ode_rk45",
"parameter vector", theta);
93 check_finite(
"integrate_ode_rk45",
"continuous data", x);
98 check_less(
"integrate_ode_rk45",
"initial time", t0_dbl, ts_dbl[0]);
100 if (relative_tolerance <= 0)
102 relative_tolerance,
"",
", must be greater than 0");
103 if (absolute_tolerance <= 0)
105 absolute_tolerance,
"",
", must be greater than 0");
106 if (max_num_steps <= 0)
108 ", must be greater than 0");
114 std::vector<double> ts_vec(ts.size() + 1);
116 std::copy(ts_dbl.begin(), ts_dbl.end(), ts_vec.begin() + 1);
118 std::vector<std::vector<typename stan::return_type<T1, T2, T_t0, T_ts>::type>>
120 coupled_ode_observer<F, T1, T2, T_t0, T_ts> observer(f, y0, theta, t0, ts, x,
122 bool observer_initial_recorded =
false;
126 auto filtered_observer
127 = [&](
const std::vector<double>& coupled_state,
double t) ->
void {
128 if (!observer_initial_recorded) {
129 observer_initial_recorded =
true;
132 observer(coupled_state, t);
136 std::vector<double> initial_coupled_state = coupled_system.initial_state();
138 const double step_size = 0.1;
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));
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.
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.
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.
const kernel_cl< in_buffer, out_buffer, int, int > copy("copy", {indexing_helpers, copy_kernel_code})
See the docs for copy() .