1 #ifndef STAN_MATH_REV_MAT_FUNCTOR_INTEGRATOR_DAE_HPP
2 #define STAN_MATH_REV_MAT_FUNCTOR_INTEGRATOR_DAE_HPP
3
4 #include <stan/math/rev/mat/functor/idas_forward_system.hpp>
5 #include <stan/math/rev/mat/functor/idas_integrator.hpp>
6 #include <ostream>
7 #include <vector>
8
9 namespace stan {
10 namespace math {
11 /**
12 * Return the solutions for a semi-explicit DAE system with residual
13 * specified by functor F,
14 * given the specified consistent initial state yy0 and yp0.
15 *
16 * @tparam DAE type of DAE system
17 * @tparam Tpar scalar type of parameter theta
18 * @param[in] f functor for the base ordinary differential equation
19 * @param[in] yy0 initial state
20 * @param[in] yp0 initial derivative state
21 * @param[in] t0 initial time
22 * @param[in] ts times of the desired solutions, in strictly
23 * increasing order, all greater than the initial time
24 * @param[in] theta parameters
25 * @param[in] x_r real data
26 * @param[in] x_i int data
27 * @param[in] rtol relative tolerance passed to IDAS, requred <10^-3
28 * @param[in] atol absolute tolerance passed to IDAS, problem-dependent
29 * @param[in] max_num_steps maximal number of admissable steps
30 * between time-points
31 * @param[in] msgs message
32 * @return a vector of states, each state being a vector of the
33 * same size as the state variable, corresponding to a time in ts.
34 */
35 template <typename F, typename Tpar>
integrate_dae(const F & f,const std::vector<double> & yy0,const std::vector<double> & yp0,double t0,const std::vector<double> & ts,const std::vector<Tpar> & theta,const std::vector<double> & x_r,const std::vector<int> & x_i,const double rtol,const double atol,const int64_t max_num_steps=idas_integrator::IDAS_MAX_STEPS,std::ostream * msgs=nullptr)36 std::vector<std::vector<Tpar> > integrate_dae(
37 const F& f, const std::vector<double>& yy0, const std::vector<double>& yp0,
38 double t0, const std::vector<double>& ts, const std::vector<Tpar>& theta,
39 const std::vector<double>& x_r, const std::vector<int>& x_i,
40 const double rtol, const double atol,
41 const int64_t max_num_steps = idas_integrator::IDAS_MAX_STEPS,
42 std::ostream* msgs = nullptr) {
43 /* it doesn't matter here what values @c eq_id has, as we
44 don't allow yy0 or yp0 to be parameters */
45 const std::vector<int> dummy_eq_id(yy0.size(), 0);
46
47 stan::math::idas_integrator solver(rtol, atol, max_num_steps);
48 stan::math::idas_forward_system<F, double, double, Tpar> dae{
49 f, dummy_eq_id, yy0, yp0, theta, x_r, x_i, msgs};
50
51 dae.check_ic_consistency(t0, atol);
52
53 return solver.integrate(dae, t0, ts);
54 }
55 } // namespace math
56 } // namespace stan
57
58 #endif
59