24 #ifndef O2SCL_ODE_BV_MULTISHOOT_H 25 #define O2SCL_ODE_BV_MULTISHOOT_H 32 #include <o2scl/astep.h> 33 #include <o2scl/astep_gsl.h> 34 #include <o2scl/ode_iv_solve.h> 35 #include <o2scl/gsl_mroot_hybrids.h> 37 #ifndef DOXYGEN_NO_O2NS 48 template <
class func_t=ode_funct<>,
49 class vec_t=ubvector,
class alloc_vec_t=ubvector,
50 class alloc_t=ubvector_alloc,
class vec_
int_t=ubvector_
int_base,
64 virtual int solve(vec_t &mesh,
int &n_func, vec_t &y_start,
65 func_t &left_b, func_t &right_b, func_t &extra_b,
66 func_t &derivs, vec_t &x_save, mat_t &y_save) {
80 int nsolve=y_start.size();
86 alloc_vec_t,alloc_t,vec_int_t> >
88 alloc_vec_t,alloc_t,vec_int_t>::
solve_fun);
92 int ret=this->
mrootp->msolve(nsolve,sx,mfm);
115 gsl_mroot_hybrids<mm_funct<> > def_mroot;
117 #ifndef DOXYGEN_INTERNAL 152 ubvector y((*this->l_n_func)),y2((*this->l_n_func));
156 for(
size_t i=0;i<(*this->
l_y_start).size();i++) {
161 for(
size_t k=0;k<(*this->
l_mesh).size()-1;k++) {
169 (*this->
l_left_b)(xa,(*this->l_n_func),sx,y);
171 for(
int i=0;i<(*this->
l_n_func);i++)
172 y[i]=sx[i+(*this->l_n_func)*(1+k)];
178 int ngrid=((*this->
l_x_save).size()-1)/((*this->l_mesh).size()-1)+1;
180 ubmatrix yysave(ngrid,(*this->l_n_func));
182 if (k!=((*this->l_mesh).size()-2)) {
183 xb=(*this->
l_mesh)[k+1]-((*this->l_mesh)[k+1]-
184 (*this->
l_mesh)[k])/ngrid;
187 this->oisp->
solve_grid(xa,xb,h,(*this->l_n_func),y,ngrid,
188 xxsave,yysave,(*this->l_derivs));
190 for(
int i=0;i<ngrid;i++) {
191 (*this->
l_x_save)[i+k*(ngrid)]=xxsave[i];
192 for(
int j=0;j<(*this->
l_n_func);j++) {
193 (*this->
l_y_save)[i+k*(ngrid)][j]=yysave[i][j];
201 (xa,xb,h,(*this->l_n_func),y,y2,(*this->l_derivs));
206 if (k==(*this->l_mesh).size()-2) {
207 (*this->
l_right_b)(xb,(*this->l_n_func),sx,y);
209 for(
int i=0;i<(*this->
l_n_func);i++) {
215 for(
int i=0;i<(*this->
l_n_func);i++) {
217 sy[i+k*(*this->
l_n_func)]= y2[i]-y[i];
223 (*this->
l_extra_b)(xb,(*this->l_n_func),sx,y);
224 for(
int i=0;i<(*this->
l_n_func);i++) {
225 sy[i+(int((*this->l_mesh).size()-1))*(*this->l_n_func)]=y[i];
235 #ifndef DOXYGEN_NO_O2NS
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
std::function< int(size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &) > mm_funct
Array of multi-dimensional functions typedef.
ode_iv_solve< func_t, vec_t, alloc_vec_t, alloc_t > * oisp
The initial value solver.
Solve a ODE boundary value problem by multishooting.
gsl_mroot_hybrids< mm_funct<> > * mrootp
The equation solver.
Multidimensional root-finding [abstract base].
One-dimensional solver [abstract base].
int solve_final_value(double x0, double x1, double h, size_t n, vec_t &ystart, vec_t ¥d, func_t &derivs)
Solve the initial-value problem to get the final value.
int solve_fun(size_t nv, const vec_t &sx, vec_t &sy)
Function to solve.
int solve_grid(double h, size_t n, size_t nsol, vec_t &xsol, mat_t &ysol, mat_t &err_sol, mat_t &dydx_sol, func_t &derivs)
Solve the initial-value problem from x0 to x1 over a grid storing derivatives and errors...