Automatic Differentiation
 
Loading...
Searching...
No Matches
ode_rk45.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_FUNCTOR_ODE_RK45_HPP
2#define STAN_MATH_PRIM_FUNCTOR_ODE_RK45_HPP
3
9#include <boost/numeric/odeint.hpp>
10#include <ostream>
11#include <tuple>
12#include <vector>
13
14namespace stan {
15namespace math {
16
55template <typename F, typename T_y0, typename T_t0, typename T_ts,
56 typename... Args, require_eigen_vector_t<T_y0>* = nullptr>
57std::vector<Eigen::Matrix<stan::return_type_t<T_y0, T_t0, T_ts, Args...>,
58 Eigen::Dynamic, 1>>
59ode_rk45_tol_impl(const char* function_name, const F& f, const T_y0& y0_arg,
60 T_t0 t0, const std::vector<T_ts>& ts,
61 double relative_tolerance, double absolute_tolerance,
62 long int max_num_steps, // NOLINT(runtime/int)
63 std::ostream* msgs, const Args&... args) {
64 using boost::numeric::odeint::integrate_times;
65 using boost::numeric::odeint::make_dense_output;
66 using boost::numeric::odeint::max_step_checker;
67 using boost::numeric::odeint::no_progress_error;
68 using boost::numeric::odeint::runge_kutta_dopri5;
69 using boost::numeric::odeint::vector_space_algebra;
70
71 using T_y0_t0 = return_type_t<T_y0, T_t0>;
72
73 Eigen::Matrix<T_y0_t0, Eigen::Dynamic, 1> y0
74 = y0_arg.template cast<T_y0_t0>();
75
76 check_finite(function_name, "initial state", y0);
77 check_finite(function_name, "initial time", t0);
78 check_finite(function_name, "times", ts);
79
80 std::tuple<ref_type_t<Args>...> args_ref_tuple(args...);
81
83 [&](const auto&... args_ref) {
84 // Code from https://stackoverflow.com/a/17340003
85 std::vector<int> unused_temp{
86 0,
87 (check_finite(function_name, "ode parameters and data", args_ref),
88 0)...};
89 },
90 args_ref_tuple);
91
92 check_nonzero_size(function_name, "initial state", y0);
93 check_nonzero_size(function_name, "times", ts);
94 check_sorted(function_name, "times", ts);
95 check_less(function_name, "initial time", t0, ts[0]);
96
97 check_positive_finite(function_name, "relative_tolerance",
98 relative_tolerance);
99 check_positive_finite(function_name, "absolute_tolerance",
100 absolute_tolerance);
101 check_positive(function_name, "max_num_steps", max_num_steps);
102
103 using return_t = return_type_t<T_y0, T_t0, T_ts, Args...>;
104 // creates basic or coupled system by template specializations
105 auto&& coupled_system = math::apply(
106 [&](const auto&... args_ref) {
108 args_ref...);
109 },
110 args_ref_tuple);
111
112 // first time in the vector must be time of initial state
113 std::vector<double> ts_vec(ts.size() + 1);
114 ts_vec[0] = value_of(t0);
115 for (size_t i = 0; i < ts.size(); ++i)
116 ts_vec[i + 1] = value_of(ts[i]);
117
118 std::vector<Eigen::Matrix<return_t, Eigen::Dynamic, 1>> y;
119 y.reserve(ts.size());
120 bool observer_initial_recorded = false;
121 size_t time_index = 0;
122
123 // avoid recording of the initial state which is included by the
124 // conventions of odeint in the output
125 auto filtered_observer
126 = [&](const std::vector<double>& coupled_state, double t) -> void {
127 if (!observer_initial_recorded) {
128 observer_initial_recorded = true;
129 return;
130 }
132 [&](const auto&... args_ref) {
133 y.emplace_back(ode_store_sensitivities(
134 f, coupled_state, y0, t0, ts[time_index], msgs, args_ref...));
135 },
136 args_ref_tuple);
137 time_index++;
138 };
139
140 // the coupled system creates the coupled initial state
141 std::vector<double> initial_coupled_state = coupled_system.initial_state();
142
143 const double step_size = 0.1;
144 try {
145 integrate_times(
146 make_dense_output(absolute_tolerance, relative_tolerance,
147 runge_kutta_dopri5<std::vector<double>, double,
148 std::vector<double>, double>()),
149 std::ref(coupled_system), initial_coupled_state, std::begin(ts_vec),
150 std::end(ts_vec), step_size, filtered_observer,
151 max_step_checker(max_num_steps));
152 } catch (const no_progress_error& e) {
153 throw_domain_error(function_name, "", ts_vec[time_index + 1],
154 "Failed to integrate to next output time (",
155 ") in less than max_num_steps steps");
156 }
157
158 return y;
159}
160
197template <typename F, typename T_y0, typename T_t0, typename T_ts,
198 typename... Args, require_eigen_vector_t<T_y0>* = nullptr>
199std::vector<Eigen::Matrix<stan::return_type_t<T_y0, T_t0, T_ts, Args...>,
200 Eigen::Dynamic, 1>>
201ode_rk45_tol(const F& f, const T_y0& y0_arg, T_t0 t0,
202 const std::vector<T_ts>& ts, double relative_tolerance,
203 double absolute_tolerance,
204 long int max_num_steps, // NOLINT(runtime/int)
205 std::ostream* msgs, const Args&... args) {
206 return ode_rk45_tol_impl("ode_rk45_tol", f, y0_arg, t0, ts,
207 relative_tolerance, absolute_tolerance,
208 max_num_steps, msgs, args...);
209}
210
243template <typename F, typename T_y0, typename T_t0, typename T_ts,
244 typename... Args, require_eigen_vector_t<T_y0>* = nullptr>
245std::vector<Eigen::Matrix<stan::return_type_t<T_y0, T_t0, T_ts, Args...>,
246 Eigen::Dynamic, 1>>
247ode_rk45(const F& f, const T_y0& y0, T_t0 t0, const std::vector<T_ts>& ts,
248 std::ostream* msgs, const Args&... args) {
249 double relative_tolerance = 1e-6;
250 double absolute_tolerance = 1e-6;
251 long int max_num_steps = 1e6; // NOLINT(runtime/int)
252
253 return ode_rk45_tol_impl("ode_rk45", f, y0, t0, ts, relative_tolerance,
254 absolute_tolerance, max_num_steps, msgs, args...);
255}
256
257} // namespace math
258} // namespace stan
259#endif
require_t< is_eigen_vector< std::decay_t< T > > > require_eigen_vector_t
Require type satisfies is_eigen_vector.
typename return_type< Ts... >::type return_type_t
Convenience type for the return type of the specified template parameters.
std::vector< Eigen::Matrix< stan::return_type_t< T_y0, T_t0, T_ts, Args... >, Eigen::Dynamic, 1 > > ode_rk45_tol_impl(const char *function_name, const F &f, const T_y0 &y0_arg, T_t0 t0, const std::vector< T_ts > &ts, double relative_tolerance, double absolute_tolerance, long int max_num_steps, std::ostream *msgs, const Args &... args)
Solve the ODE initial value problem y' = f(t, y), y(t0) = y0 at a set of times, { t1,...
Definition ode_rk45.hpp:59
static constexpr double e()
Return the base of the natural logarithm.
Definition constants.hpp:20
std::vector< Eigen::Matrix< stan::return_type_t< T_y0, T_t0, T_ts, Args... >, Eigen::Dynamic, 1 > > ode_rk45(const F &f, const T_y0 &y0, T_t0 t0, const std::vector< T_ts > &ts, std::ostream *msgs, const Args &... args)
Solve the ODE initial value problem y' = f(t, y), y(t0) = y0 at a set of times, { t1,...
Definition ode_rk45.hpp:247
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition value_of.hpp:18
void throw_domain_error(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
Throw a domain error with a consistently formatted message.
std::vector< Eigen::Matrix< stan::return_type_t< T_y0, T_t0, T_ts, Args... >, Eigen::Dynamic, 1 > > ode_rk45_tol(const F &f, const T_y0 &y0_arg, T_t0 t0, const std::vector< T_ts > &ts, double relative_tolerance, double absolute_tolerance, long int max_num_steps, std::ostream *msgs, const Args &... args)
Solve the ODE initial value problem y' = f(t, y), y(t0) = y0 at a set of times, { t1,...
Definition ode_rk45.hpp:201
void check_finite(const char *function, const char *name, const T_y &y)
Return true if all values in y are finite.
Eigen::VectorXd ode_store_sensitivities(const F &f, const std::vector< double > &coupled_state, const Eigen::Matrix< T_y0_t0, Eigen::Dynamic, 1 > &y0, T_t0 t0, T_t t, std::ostream *msgs, const Args &... args)
When all arguments are arithmetic, there are no sensitivities to store, so the function just returns ...
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_positive(const char *function, const char *name, const T_y &y)
Check if y is positive.
void check_sorted(const char *function, const char *name, const EigVec &y)
Check if the specified vector is sorted into increasing order (repeated values are okay).
void check_less(const char *function, const char *name, const T_y &y, const T_high &high, Idxs... idxs)
Throw an exception if y is not strictly less than high.
constexpr decltype(auto) apply(F &&f, Tuple &&t, PreArgs &&... pre_args)
Definition apply.hpp:52
void check_positive_finite(const char *function, const char *name, const T_y &y)
Check if y is positive and finite.
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...