include/continewz/ode_series.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : ode_series.hpp
* DESCRIPTION: Series solutions to dynamical systems
* COPYRIGHT  : (C) 2011  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX_ODE_SERIES_HPP
#define __MMX_ODE_SERIES_HPP
#include <algebramix/series_vector.hpp>
#include <algebramix/series_matrix.hpp>
#include <multimix/slp_multivariate.hpp>

namespace mmx {
#define Homotopy_variant(C) typename homotopy_variant_helper<C>::HV
#define TMPL template<typename C>
#define Series series<C>
#define Series_rep series_rep<C>
#define Function_C slp_tangent<C>
#define Function_Series slp_tangent<Series >

/******************************************************************************
* Integration of dynamical systems
******************************************************************************/

template<typename C>
class solve_ode_series_rep: public recursive_series_rep<vector<C> > {
  Function_Series fun;
  vector<C>       c;
public:
  solve_ode_series_rep (const Function_Series& fun2, const vector<C>& c2):
    recursive_series_rep<vector<C> > (CF(c2)), fun (fun2), c (c2) {}
  syntactic expression (const syntactic& z) const {
    (void) z;
    return apply ("solve_ode", flatten (fun), flatten (c)); }
  series<vector<C> > initialize () {
    this->initial (0)= c;
    vector<series<C> > in = as_vector (this->me ());
    vector<series<C> > out= integrate (eval (fun, in));
    return from_vector (out); }
};

template<typename C> inline vector<Series >
solve_ode (const Function_Series& fun, const vector<C>& c) {
  series_rep<vector<C> >* rep= new solve_ode_series_rep<C> (fun, c);
  return as_vector (recursive (series<vector<C> > (rep)));
}

template<typename C> inline vector<Series >
solve_ode_series (const Function_C& fun, const vector<C>& c) {
  return solve_ode (as<Function_Series > (fun), c);
}

template<typename C> vector<generic>
solve_ode_series (const vector<multivariate<sparse_polynomial<rational> > >& f,
                  const vector<multivariate_coordinate<> >& vars,
                  const vector<C>& c) {
  Function_C fun= as_slp_tangent (f, vars, CF (c));
  return as<vector<generic> > (solve_ode_series (fun, c));
}

/******************************************************************************
* First variation at solution
******************************************************************************/

template<typename C> inline matrix<Series >
first_variation (const Function_Series& fun, const vector<Series >& at) {
  matrix<Series > jac= jacobian (fun, at);
  return solve_lde (jac);
}

template<typename C> inline matrix<Series >
first_variation_series (const Function_C& fun, const vector<Series >& at) {
  return first_variation (as<Function_Series > (fun), at);
}

template<typename C> inline matrix<Series >
first_variation_series (const Function_C& fun, const vector<C>& c) {
  vector<Series > at= solve_ode_series (fun, c);
  return first_variation_series (fun, at);
}

template<typename C> matrix<generic>
first_variation_series (const vector<multivariate<sparse_polynomial<rational> > >& f,
                        const vector<multivariate_coordinate<> >& vars,
                        const vector<C>& c) {
  Function_C fun= as_slp_tangent (f, vars, CF (c));
  return as<matrix<generic> > (first_variation_series (fun, c));
}

} // namespace mmx
#undef TMPL
#undef Series
#undef Series_rep
#undef Function_C
#undef Function_Series
#endif // __MMX_ODE_SERIES_HPP

---