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 <ostream>
19 #include <vector>
20 
21 namespace stan {
22 namespace math {
23 
24 /**
25  * Return the solutions for the specified system of ordinary
26  * differential equations given the specified initial state,
27  * initial times, times of desired solution, and parameters and
28  * data, writing error and warning messages to the specified
29  * stream.
30  *
31  * <b>Warning:</b> If the system of equations is stiff, roughly
32  * defined by having varying time scales across dimensions, then
33  * this solver is likely to be slow.
34  *
35  * This function is templated to allow the initial times to be
36  * either data or autodiff variables and the parameters to be data
37  * or autodiff variables. The autodiff-based implementation for
38  * reverse-mode are defined in namespace <code>stan::math</code>
39  * and may be invoked via argument-dependent lookup by including
40  * their headers.
41  *
42  * This function uses the <a
43  * href="http://en.wikipedia.org/wiki/Dormand–Prince_method">Dormand-Prince
44  * method</a> as implemented in Boost's <code>
45  * boost::numeric::odeint::runge_kutta_dopri5</code> integrator.
46  *
47  * @tparam F type of ODE system function.
48  * @tparam T1 type of scalars for initial values.
49  * @tparam T2 type of scalars for parameters.
50  * @param[in] f functor for the base ordinary differential equation.
51  * @param[in] y0 initial state.
52  * @param[in] t0 initial time.
53  * @param[in] ts times of the desired solutions, in strictly
54  * increasing order, all greater than the initial time.
55  * @param[in] theta parameter vector for the ODE.
56  * @param[in] x continuous data vector for the ODE.
57  * @param[in] x_int integer data vector for the ODE.
58  * @param[out] msgs the print stream for warning messages.
59  * @param[in] relative_tolerance relative tolerance parameter
60  * for Boost's ode solver. Defaults to 1e-6.
61  * @param[in] absolute_tolerance absolute tolerance parameter
62  * for Boost's ode solver. Defaults to 1e-6.
63  * @param[in] max_num_steps maximum number of steps to take within
64  * the Boost ode solver.
65  * @return a vector of states, each state being a vector of the
66  * same size as the state variable, corresponding to a time in ts.
67  */
68 template <typename F, typename T1, typename T2>
70 integrate_ode_rk45(const F& f, const std::vector<T1>& y0, double t0,
71  const std::vector<double>& ts, const std::vector<T2>& theta,
72  const std::vector<double>& x, const std::vector<int>& x_int,
73  std::ostream* msgs = nullptr,
74  double relative_tolerance = 1e-6,
75  double absolute_tolerance = 1e-6, int max_num_steps = 1E6) {
76  using boost::numeric::odeint::integrate_times;
77  using boost::numeric::odeint::make_dense_output;
78  using boost::numeric::odeint::max_step_checker;
79  using boost::numeric::odeint::runge_kutta_dopri5;
80 
81  check_finite("integrate_ode_rk45", "initial state", y0);
82  check_finite("integrate_ode_rk45", "initial time", t0);
83  check_finite("integrate_ode_rk45", "times", ts);
84  check_finite("integrate_ode_rk45", "parameter vector", theta);
85  check_finite("integrate_ode_rk45", "continuous data", x);
86 
87  check_nonzero_size("integrate_ode_rk45", "times", ts);
88  check_nonzero_size("integrate_ode_rk45", "initial state", y0);
89  check_ordered("integrate_ode_rk45", "times", ts);
90  check_less("integrate_ode_rk45", "initial time", t0, ts[0]);
91 
92  if (relative_tolerance <= 0)
93  invalid_argument("integrate_ode_rk45", "relative_tolerance,",
94  relative_tolerance, "", ", must be greater than 0");
95  if (absolute_tolerance <= 0)
96  invalid_argument("integrate_ode_rk45", "absolute_tolerance,",
97  absolute_tolerance, "", ", must be greater than 0");
98  if (max_num_steps <= 0)
99  invalid_argument("integrate_ode_rk45", "max_num_steps,", max_num_steps, "",
100  ", must be greater than 0");
101 
102  // creates basic or coupled system by template specializations
103  coupled_ode_system<F, T1, T2> coupled_system(f, y0, theta, x, x_int, msgs);
104 
105  // first time in the vector must be time of initial state
106  std::vector<double> ts_vec(ts.size() + 1);
107  ts_vec[0] = t0;
108  for (size_t n = 0; n < ts.size(); n++)
109  ts_vec[n + 1] = ts[n];
110 
111  std::vector<std::vector<double> > y_coupled(ts_vec.size());
112  coupled_ode_observer observer(y_coupled);
113 
114  // the coupled system creates the coupled initial state
115  std::vector<double> initial_coupled_state = coupled_system.initial_state();
116 
117  const double step_size = 0.1;
118  integrate_times(
119  make_dense_output(absolute_tolerance, relative_tolerance,
120  runge_kutta_dopri5<std::vector<double>, double,
121  std::vector<double>, double>()),
122  boost::ref(coupled_system), initial_coupled_state, boost::begin(ts_vec),
123  boost::end(ts_vec), step_size, observer, max_step_checker(max_num_steps));
124 
125  // remove the first state corresponding to the initial value
126  y_coupled.erase(y_coupled.begin());
127 
128  // the coupled system also encapsulates the decoupling operation
129  return coupled_system.decouple_states(y_coupled);
130 }
131 
132 } // namespace math
133 
134 } // namespace stan
135 
136 #endif
#define F(x, y, z)
void check_finite(const char *function, const char *name, const T_y &y)
void check_nonzero_size(const char *function, const char *name, const T_y &y)
void check_ordered(const char *function, const char *name, const std::vector< T_y > &y)
std::vector< std::vector< typename stan::return_type< T1, T2 >::type > > integrate_ode_rk45(const F &f, const std::vector< T1 > &y0, double t0, const std::vector< double > &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)
void invalid_argument(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
double e()
Definition: constants.hpp:89
void check_less(const char *function, const char *name, const T_y &y, const T_high &high)
Definition: check_less.hpp:72