24 #ifndef O2SCL_ODE_BV_MSHOOT_H 25 #define O2SCL_ODE_BV_MSHOOT_H 31 #include <o2scl/ode_bv_solve.h> 33 #ifndef DOXYGEN_NO_O2NS 72 vec_t &x_bound, mat_t &y_bound,
73 vec_int_t &index, func_t &derivs) {
84 int lhs_unks=0, rhs_conts=0;
88 for (
size_t i=0;i<n;i++) {
89 if (index[i]<2) lhs_unks++;
90 if (index[i]%2==1) rhs_conts++;
94 if (lhs_unks!=rhs_conts) {
95 O2SCL_ERR2(
"Incorrect boundary conditions in ",
96 "ode_bv_mshoot::solve()",gsl_einval);
100 int n_solve=lhs_unks+n*(n_bound-2);
107 for (
size_t i=0;i<n;i++) {
113 for(
size_t k=1;k<n_bound-1;k++) {
114 for(
size_t i=0;i<n;i++) {
131 mm_funct_mfptr<ode_bv_mshoot<func_t,mat_t,vec_t,vec_int_t> >
137 O2SCL_ERR(
"Solver failed in ode_bv_mshoot::solve().",ret);
146 template<
class mat_row_t>
147 int solve_store(
double h,
size_t n,
size_t n_bound, vec_t &x_bound,
148 mat_t &y_bound, vec_int_t &index,
size_t &n_sol,
149 vec_t &x_sol, mat_t &y_sol, mat_t &dydx_sol,
150 mat_t &yerr_sol, func_t &derivs) {
153 O2SCL_ERR2(
"Not enough boundaries (must be at least two) in ",
154 "ode_bv_mshoot::solve_store().",gsl_einval);
157 O2SCL_ERR2(
"Not enough room to store boundaries in ",
158 "ode_bv_mshoot::solve_store().",gsl_einval);
165 ubvector_int inxs(n_bound);
166 for(
size_t k=0;k<n_bound;k++) {
167 inxs[k]=((size_t)(((
double)k)/((double)n_bound)*
168 (((double)(n_sol))-1.0+1.0e-12)));
169 std::cout << k <<
" " << inxs[k] <<
" " << n_sol << std::endl;
172 for(
size_t k=1;k<n_bound-1;k++) {
173 if (inxs[k]==inxs[k-1] || inxs[k]==inxs[k+1]) {
174 O2SCL_ERR2(
"Not enough room to store boundaries in ",
175 "ode_bv_mshoot::solve_store().",gsl_einval);
180 for(
size_t k=0;k<n_bound-1;k++) {
181 size_t n_sol_tmp=inxs[k+1];
182 mat_row_t ystart(y_bound,k);
184 std::cout <<
"Old boundaries: " << inxs << std::endl;
186 oisp->template solve_store<mat_t,mat_row_t>
187 (x_bound[k],x_bound[k+1],h,n,ystart,
188 n_sol_tmp,x_sol,y_sol,dydx_sol,yerr_sol,derivs,inxs[k]);
190 std::cout <<
"New boundaries: " << n_sol_tmp << std::endl;
194 if (((
int)n_sol_tmp)<inxs[k+1]) {
195 for(
size_t k2=k+1;k2<n_bound;k2++) {
196 inxs[k2]=((size_t)(((
double)k2-k-1)/((
double)(n_bound-k-1))*
197 (((double)(n_sol-n_sol_tmp)))))+n_sol_tmp-1;
201 std::cout <<
"New boundaries: " << inxs << std::endl;
204 n_sol=inxs[n_bound-1];
228 #ifndef DOXYGEN_INTERNAL 264 vec_t sy, sy2, syerr, sdydx, sdydx2;
281 for(
size_t k=0;k<l_nbound-1;k++) {
283 for(
size_t i=0;i<this->
l_n;i++) {
284 if ((*this->l_index)[i]<2) {
290 for(
size_t i=0;i<this->
l_n;i++) {
298 for(
size_t k=0;k<l_nbound-1;k++) {
307 for(
size_t i=0;i<this->
l_n;i++) {
315 sy2,*this->l_derivs);
319 for(
size_t i=0;i<this->
l_n;i++) {
320 ty[j_y]=sy2[i]-yb[k+1][i];
325 for(
size_t i=0;i<this->
l_n;i++) {
326 if ((*this->l_index)[i]%2==1) {
327 double yright=yb[k+1][i];
331 ty[j_y]=(sy2[i]-yright)/yright;
348 #ifndef DOXYGEN_NO_O2NS Solve an initial-value ODE problems given an adaptive ODE stepper.
size_t mem_size
Size of recent allocation.
mat_t * l_ybound
Storage for the ending vector.
double l_h
Storage for the stepsize.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
vec_t * l_xbound
Storage for the starting vector.
int set_mroot(mroot< mm_funct<> > &root)
Set the equation solver.
ode_iv_solve< func_t, vec_t > * oisp
The solver for the initial value problem.
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 > def_ois
The default initial value solver.
virtual int msolve(size_t n, vec_t &x, func_t &func)=0
Solve func using x as an initial guess, returning x.
int solve_store(double h, size_t n, size_t n_bound, vec_t &x_bound, mat_t &y_bound, vec_int_t &index, size_t &n_sol, vec_t &x_sol, mat_t &y_sol, mat_t &dydx_sol, mat_t &yerr_sol, func_t &derivs)
Solve the boundary-value problem and store the solution.
mroot< mm_funct<> > * mrootp
The equation solver.
Multidimensional root-finding [abstract base].
Solve boundary-value ODE problems by multishooting with a generic nonlinear solver.
size_t l_nbound
The number of boundaries.
func_t * l_derivs
The functions to integrate.
One-dimensional solver [abstract base].
vec_int_t * l_index
The index defining the boundary conditions.
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.
std::function< int(double, size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &)> ode_funct
Ordinary differential equation function.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
int solve_final_value(double h, size_t n, size_t n_bound, vec_t &x_bound, mat_t &y_bound, vec_int_t &index, func_t &derivs)
Solve the boundary-value problem and store the solution.
size_t l_lhs_unks
The number of unknowns on the LHS.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
gsl_mroot_hybrids< mm_funct<> > def_mroot
The default equation solver.
size_t l_n
The number of functions.
int solve_fun(size_t nv, const vec_t &tx, vec_t &ty)
The shooting function to be solved by the multidimensional solver.
int set_iv(ode_iv_solve< func_t, vec_t > &ois)
Set initial value solver.
static const double x1[5]
Base class for boundary-value ODE solvers.