| Stan Math Library 2.11.0 reverse mode automatic differentiation

- Main Page
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

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
4 #include <stan/math/prim/arr/err/check_nonzero_size.hpp>
5 #include <stan/math/prim/arr/err/check_ordered.hpp>
6 #include <stan/math/prim/arr/functor/coupled_ode_system.hpp>
7 #include <stan/math/prim/arr/functor/coupled_ode_observer.hpp>
8 #include <stan/math/prim/scal/fun/value_of.hpp>
9 #include <stan/math/prim/scal/err/check_less.hpp>
10 #include <stan/math/prim/scal/err/check_finite.hpp>
11 #include <stan/math/prim/scal/err/invalid_argument.hpp>
12 #include <stan/math/prim/scal/meta/return_type.hpp>
13 #include <boost/numeric/odeint.hpp>
14 #include <ostream>
15 #include <vector>
16
17 namespace stan {
18
19 namespace math {
20
65 template <typename F, typename T1, typename T2>
66 std::vector<std::vector<typename stan::return_type<T1, T2>::type> >
67 integrate_ode_rk45(const F& f,
68 const std::vector<T1> y0,
69 const double t0,
70 const std::vector<double>& ts,
71 const std::vector<T2>& theta,
72 const std::vector<double>& x,
73 const std::vector<int>& x_int,
74 std::ostream* msgs = 0,
75 double relative_tolerance = 1e-6,
76 double absolute_tolerance = 1e-6,
77 int max_num_steps = 1E6) {
78 using boost::numeric::odeint::integrate_times;
79 using boost::numeric::odeint::make_dense_output;
80 using boost::numeric::odeint::runge_kutta_dopri5;
81 using boost::numeric::odeint::max_step_checker;
82
83 check_finite("integrate_ode_rk45", "initial state", y0);
84 check_finite("integrate_ode_rk45", "initial time", t0);
85 check_finite("integrate_ode_rk45", "times", ts);
86 check_finite("integrate_ode_rk45", "parameter vector", theta);
87 check_finite("integrate_ode_rk45", "continuous data", x);
88
89 check_nonzero_size("integrate_ode_rk45", "times", ts);
90 check_nonzero_size("integrate_ode_rk45", "initial state", y0);
91 check_ordered("integrate_ode_rk45", "times", ts);
92 check_less("integrate_ode_rk45", "initial time", t0, ts[0]);
93
94 if (relative_tolerance <= 0)
95 invalid_argument("integrate_ode_rk45",
96 "relative_tolerance,", relative_tolerance,
97 "", ", must be greater than 0");
98 if (absolute_tolerance <= 0)
99 invalid_argument("integrate_ode_rk45",
100 "absolute_tolerance,", absolute_tolerance,
101 "", ", must be greater than 0");
102 if (max_num_steps <= 0)
103 invalid_argument("integrate_ode_rk45",
104 "max_num_steps,", max_num_steps,
105 "", ", must be greater than 0");
106
107 // creates basic or coupled system by template specializations
108 coupled_ode_system<F, T1, T2>
109 coupled_system(f, y0, theta, x, x_int, msgs);
110
111 // first time in the vector must be time of initial state
112 std::vector<double> ts_vec(ts.size() + 1);
113 ts_vec[0] = t0;
114 for (size_t n = 0; n < ts.size(); n++)
115 ts_vec[n+1] = ts[n];
116
117 std::vector<std::vector<double> > y_coupled(ts_vec.size());
118 coupled_ode_observer observer(y_coupled);
119
120 // the coupled system creates the coupled initial state
121 std::vector<double> initial_coupled_state
122 = coupled_system.initial_state();
123
124 const double step_size = 0.1;
125 integrate_times(make_dense_output(absolute_tolerance,
126 relative_tolerance,
127 runge_kutta_dopri5<std::vector<double>,
128 double,
129 std::vector<double>,
130 double>() ),
131 coupled_system,
132 initial_coupled_state,
133 boost::begin(ts_vec), boost::end(ts_vec),
134 step_size,
135 observer,
136 max_step_checker(max_num_steps));
137
138 // remove the first state corresponding to the initial value
139 y_coupled.erase(y_coupled.begin());
140
141 // the coupled system also encapsulates the decoupling operation
142 return coupled_system.decouple_states(y_coupled);
143 }
144
145 }
146
147 }
148
149 #endif
stan::math::check_less
bool check_less(const char *function, const char *name, const T_y &y, const T_high &high)
Return true if y is strictly less than high.
Definition: check_less.hpp:81
stan
Definition: log_sum_exp.hpp:8
invalid_argument.hpp
return_type.hpp
coupled_ode_observer.hpp
stan::math::check_nonzero_size
bool check_nonzero_size(const char *function, const char *name, const T_y &y)
Return true if the specified matrix/vector is of non-zero size.
Definition: check_nonzero_size.hpp:30
stan::math::check_ordered
bool check_ordered(const char *function, const char *name, const std::vector< T_y > &y)
Return true if the specified vector is sorted into strictly increasing order.
Definition: check_ordered.hpp:30
stan::math::coupled_ode_observer
Observer for the coupled states.
Definition: coupled_ode_observer.hpp:15
stan::math::integrate_ode_rk45
std::vector< std::vector< typename stan::return_type< T1, T2 >::type > > integrate_ode_rk45(const F &f, const std::vector< T1 > y0, const 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=0, 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 ...
Definition: integrate_ode_rk45.hpp:67
stan::math::invalid_argument
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.
Definition: invalid_argument.hpp:31
stan::math::e
double e()
Return the base of the natural logarithm.
Definition: constants.hpp:95
stan::math::check_finite
bool check_finite(const char *function, const char *name, const T_y &y)
Return true if y is finite.
Definition: check_finite.hpp:62
check_ordered.hpp
check_finite.hpp
stan::math::coupled_ode_system
Base template class for a coupled ordinary differential equation system, which adds sensitivities to ...
Definition: coupled_ode_system.hpp:25
coupled_ode_system.hpp
value_of.hpp
check_nonzero_size.hpp
check_less.hpp

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