ode_bv_solve.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 #ifndef O2SCL_ODE_BV_SOLVE_H
24 #define O2SCL_ODE_BV_SOLVE_H
25 
26 /** \file ode_bv_solve.h
27  \brief File defining \ref o2scl::ode_bv_solve
28 */
29 
30 #include <string>
31 #include <o2scl/astep.h>
32 #include <o2scl/astep_gsl.h>
33 #include <o2scl/ode_iv_solve.h>
34 #include <o2scl/mroot_hybrids.h>
35 
36 #ifndef DOXYGEN_NO_O2NS
37 namespace o2scl {
38 #endif
39 
40  /** \brief Base class for boundary-value ODE solvers
41 
42  This class is experimental.
43  */
44  class ode_bv_solve {
45 
46  public:
47 
48  ode_bv_solve() {
49  verbose=0;
50  }
51 
52  virtual ~ode_bv_solve() {}
53 
54  /** \name Values for index arrays */
55  //@{
56  /// Unknown on both the left and right boundaries
57  static const int unk=0;
58  /// Known on the right boundary
59  static const int right=1;
60  /// Known on the left boundary
61  static const int left=2;
62  /// Known on both the left and right boundaries
63  static const int both=3;
64  //@}
65 
66  /// Set output level
67  int verbose;
68 
69  };
70 
71  /** \brief Solve boundary-value ODE problems by shooting from one
72  boundary to the other
73 
74  This class is experimental.
75 
76  Documentation links for default template arguments
77  - \c func_t - \ref ode_funct11
78  - \c vec_t - \ref boost::numeric::ublas::vector < double >
79  - \c vec_int_t - \ref boost::numeric::ublas::vector < int >
80  */
81  template<class func_t=ode_funct11,
83  class vec_int_t=boost::numeric::ublas::vector<int> >
84  class ode_bv_shoot : public ode_bv_solve {
85 
86  public:
87 
90 
91  ode_bv_shoot() {
92  oisp=&def_ois;
93  mrootp=&def_mroot;
94  mem_size=0;
95  }
96 
97  virtual ~ode_bv_shoot() {
98  }
99 
100  /// Allocate internal storage
101  void allocate(size_t n) {
102  if (n!=mem_size) {
103  sy.resize(n);
104  sy2.resize(n);
105  syerr.resize(n);
106  sdydx.resize(n);
107  mem_size=n;
108  }
109  return;
110  }
111 
112  /** \brief Solve the boundary-value problem and store the
113  solution
114 
115  Given the \c n initial values of the functions in \c ystart,
116  this function integrates the ODEs specified in \c derivs over
117  the interval from \c x0 to \c x1 with an initial stepsize of
118  \c h. The final values of the function are given in \c yend,
119  the derivatives in \c dydx_end, and the associated errors are
120  given in \c yerr. The initial values of \c yend and \c yerr
121  are ignored.
122 
123  */
124  int solve_final_value(double x0, double x1, double h, size_t n,
125  vec_t &ystart, vec_t &yend, vec_int_t &index,
126  vec_t &yerr, vec_t &dydx_end, func_t &derivs) {
127 
128  // Get pointers to inputs for access by solve function
129  this->l_index=&index;
130  this->l_ystart=&ystart;
131  this->l_yend=&yend;
132  this->l_yerr=&yerr;
133  this->l_dydx_end=&dydx_end;
134  this->l_x0=x0;
135  this->l_x1=x1;
136  this->l_h=h;
137  this->l_derivs=&derivs;
138  this->l_n=n;
139 
140  // lhs_unks is the number of unknowns on the LHS
141  // rhs_conts is the number of constraints on the RHS
142  int lhs_unks=0, rhs_conts=0;
143 
144  // Count the number of variables we'll need to solve for
145  for (size_t i=0;i<n;i++) {
146  if (index[i]<2) lhs_unks++;
147  if (index[i]%2==1) rhs_conts++;
148  }
149 
150  // Make sure that the boundary conditions make sense
151  if (lhs_unks!=rhs_conts) {
152  O2SCL_ERR2("Incorrect boundary conditions in ",
153  "ode_bv_shoot::solve()",gsl_einval);
154  }
155 
156  // Make space for the solution
157  ubvector tx(lhs_unks);
158 
159  // Copy initial guesses from ystart
160  int j=0;
161  for (size_t i=0;i<n;i++) {
162  if (index[i]<2) {
163  tx[j]=ystart[i];
164  j++;
165  }
166  }
167 
168  allocate(n);
169 
170  // The function object
171  mm_funct_mfptr<ode_bv_shoot<func_t,vec_t,vec_int_t> >
173 
174  // Solve
175  int ret=this->mrootp->msolve(lhs_unks,tx,mfm);
176  if (ret!=0) {
177  O2SCL_ERR("Solver failed in ode_bv_shoot::solve().",ret);
178  }
179 
180  // Copy the solution back to ystart
181  j=0;
182  for (size_t i=0;i<n;i++) {
183  if (index[i]<2) {
184  ystart[i]=tx[j];
185  j++;
186  }
187  }
188 
189  return 0;
190  }
191 
192  /** \brief Solve the boundary-value problem and store the
193  solution
194  */
195  template<class mat_t, class mat_row_t>
196  int solve_store(double x0, double x1, double h, size_t n,
197  vec_t &ystart, vec_t &yend,
198  vec_int_t &index, size_t &n_sol, vec_t &x_sol,
199  mat_t &y_sol, mat_t &yerr_sol, mat_t &dydx_sol,
200  func_t &derivs) {
201 
202  mat_row_t yerr(yerr_sol,0);
203  mat_row_t dydx_end(dydx_sol,0);
204 
205  // Solve for the final value
206  solve_final_value(x0,x1,h,n,ystart,yend,index,yerr,dydx_end,derivs);
207 
208  // Copy ystart to y_sol[0]
209  for(size_t i=0;i<n;i++) {
210  y_sol[0][i]=ystart[i];
211  }
212 
213  // Evaluate the function one more time to create the table
214  oisp->solve_store<mat_t,mat_row_t>(this->l_x0,this->l_x1,this->l_h,
215  this->l_n,n_sol,x_sol,
216  y_sol,yerr_sol,dydx_sol,derivs);
217 
218  // Copy the stored solution back to ystart and yend
219  for (size_t i=0;i<n;i++) {
220  if (index[i]%2==0) {
221  yend[i]=y_sol[n_sol-1][i];
222  }
223  ystart[i]=y_sol[0][i];
224  }
225 
226  return 0;
227  }
228 
229  /** \brief Set initial value solver
230  */
232  oisp=&ois;
233  return 0;
234  }
235 
236  /** \brief Set the equation solver */
237  int set_mroot(mroot<mm_funct<> > &root) {
238  mrootp=&root;
239  return 0;
240  }
241 
242  /// The default initial value solver
244 
245  /// The default equation solver
246  gsl_mroot_hybrids<mm_funct<> > def_mroot;
247 
248 #ifndef DOXYGEN_INTERNAL
249 
250  protected:
251 
252  /// The solver for the initial value problem
254 
255  /// The equation solver
257 
258  /// The index defining the boundary conditions
259  vec_int_t *l_index;
260 
261  /// Storage for the starting vector
262  vec_t *l_ystart;
263 
264  /// Storage for the ending vector
265  vec_t *l_yend;
266 
267  /// Storage for the starting vector
268  vec_t *l_yerr;
269 
270  /// Storage for the ending vector
271  vec_t *l_dydx_end;
272 
273  /// Storage for the starting point
274  double l_x0;
275 
276  /// Storage for the ending abcissa
277  double l_x1;
278 
279  /// Storage for the stepsize
280  double l_h;
281 
282  /// The functions to integrate
283  func_t *l_derivs;
284 
285  /// The number of functions
286  size_t l_n;
287 
288  /// \name Temporary storage for \ref solve_fun()
289  //@{
290  vec_t sy, sy2, syerr, sdydx;
291  //@}
292 
293  /// Size of recent allocation
294  size_t mem_size;
295 
296  /// The shooting function to be solved by the multidimensional solver
297  int solve_fun(size_t nv, const vec_t &tx, vec_t &ty) {
298 
299  int j;
300 
301  // Create leftmost point from combination of
302  // starting array and proposed solution
303  j=0;
304  for(size_t i=0;i<this->l_n;i++) {
305  if ((*this->l_index)[i]<2) {
306  sy[i]=tx[j];
307  j++;
308  } else {
309  sy[i]=(*this->l_ystart)[i];
310  }
311  }
312 
313  // Shoot across. We specify memory for the derivatives and
314  // errors to prevent ode_iv_solve from having to continuously
315  // allocate and reallocate it.
316  oisp->solve_final_value(this->l_x0,this->l_x1,this->l_h,
317  this->l_n,sy,sy2,*this->l_yerr,
318  *this->l_dydx_end,*this->l_derivs);
319 
320  j=0;
321  for(size_t i=0;i<this->l_n;i++) {
322  if ((*this->l_index)[i]%2==1) {
323  // Construct the equations from the rightmost point
324  if ((*this->l_yend)[i]==0.0) {
325  ty[j]=sy2[i];
326  } else {
327  ty[j]=(sy2[i]-(*this->l_yend)[i])/(*this->l_yend)[i];
328  }
329  j++;
330  } else {
331  // Otherwise copy the final values from y2 to *l_yend
332  (*this->l_yend)[i]=sy2[i];
333  }
334  }
335 
336  return 0;
337  }
338 
339 #endif
340 
341  };
342 
343  /** \brief Solve boundary-value ODE problems by shooting from one
344  boundary to the other on a grid
345 
346  This class is experimental.
347 
348  Default template arguments
349  - \c func_t - \ref ode_funct11
350  - \c vec_t - \ref boost::numeric::ublas::vector < double >
351  - \c vec_int_t - \ref boost::numeric::ublas::vector < int >
352  */
353  template<class mat_t=boost::numeric::ublas::matrix<double>,
354  class mat_row_t=boost::numeric::ublas::matrix_row<
355  boost::numeric::ublas::matrix<double> >,
356  class func_t=ode_funct<>,
357  class vec_t=boost::numeric::ublas::vector<double>,
358  class vec_int_t=boost::numeric::ublas::vector<int> >
360  public ode_bv_shoot<func_t,vec_t,vec_int_t> {
361 
362  public:
363 
367 
369  }
370 
371  virtual ~ode_bv_shoot_grid() {
372  }
373 
374  /// Desc
375  int solve_grid(double x0, double x1, double h, size_t n,
376  vec_t &ystart, vec_t &yend, vec_int_t &index,
377  size_t nsol, vec_t &xsol, mat_t &ysol,
378  mat_t &err_sol, mat_t &dydx_sol, func_t &derivs) {
379 
380  // lhs_unks is the number of unknowns on the LHS
381  // rhs_conts is the number of constraints on the RHS
382  int lhs_unks=0, rhs_conts=0;
383 
384  // Count the number of variables we'll need to solve for
385  for (size_t i=0;i<n;i++) {
386  if (index[i]<2) lhs_unks++;
387  if (index[i]%2==1) rhs_conts++;
388  }
389 
390  // Make sure that the boundary conditions make sense
391  if (lhs_unks!=rhs_conts) {
392  O2SCL_ERR2("Incorrect boundary conditions in ",
393  "ode_bv_shoot_grid::solve_grid()",gsl_einval);
394  }
395 
396  // Get pointers to inputs for access by solve function
397  this->l_n=n;
398  this->l_index=&index;
399  this->l_derivs=&derivs;
400  this->l_ystart=&ystart;
401  this->l_yend=&yend;
402  this->l_h=h;
403  l_nsol=nsol;
404  l_xsol=&xsol;
405  l_ysol=&ysol;
406  l_dydxsol=&dydx_sol;
407  l_errsol=&err_sol;
408 
409  // Make space for the solution
410  ubvector tx(lhs_unks);
411 
412  // Copy initial guesses from ysol
413  int j=0;
414  for (size_t i=0;i<n;i++) {
415  if (index[i]<2) {
416  tx[j]=ystart[i];
417  j++;
418  }
419  }
420 
421  // Allocate internal storage
422  this->allocate(n);
423 
424  // The function object
425  mm_funct_mfptr<ode_bv_shoot_grid<mat_t,mat_row_t,
426  func_t,vec_t,vec_int_t> >
427  mfm(this,&ode_bv_shoot_grid<mat_t,mat_row_t,func_t,vec_t,
428  vec_int_t>::solve_grid_fun);
429 
430  // Solve
431  int ret=this->mrootp->msolve(lhs_unks,tx,mfm);
432  if (ret!=0) {
433  O2SCL_ERR("Solver failed in ode_bv_shoot_grid::solve_grid().",ret);
434  }
435 
436  // Copy the solution back to ysol
437  j=0;
438  for (size_t i=0;i<n;i++) {
439  if (index[i]<2) {
440  ysol[0][i]=tx[j];
441  j++;
442  }
443  }
444 
445  return 0;
446  }
447 
448 #ifndef DOXYGEN_INTERNAL
449 
450  protected:
451 
452  /// Desc
453  size_t l_nsol;
454 
455  /// Desc
456  vec_t *l_xsol;
457 
458  /// Desc
459  mat_t *l_ysol;
460 
461  /// Desc
462  mat_t *l_dydxsol;
463 
464  /// Desc
465  mat_t *l_errsol;
466 
467  /// The shooting function to be solved by the multidimensional solver
468  int solve_grid_fun(size_t nv, const vec_t &tx, vec_t &ty) {
469 
470  int j;
471 
472  // Create leftmost point from combination of
473  // starting array and proposed solution
474  j=0;
475  for(size_t i=0;i<this->l_n;i++) {
476  if ((*this->l_index)[i]<2) {
477  (*l_ysol)[0][i]=tx[j];
478  j++;
479  }
480  }
481 
482  // Perform the solution
483  this->oisp->template solve_grid<mat_t,mat_row_t>
484  (this->l_h,this->l_n,l_nsol,*l_xsol,*l_ysol,*l_errsol,
485  *l_dydxsol,*this->l_derivs);
486 
487  // The last vector
488  mat_row_t yend2(*l_ysol,l_nsol-1);
489 
490  j=0;
491  for(size_t i=0;i<this->l_n;i++) {
492  if ((*this->l_index)[i]%2==1) {
493  // Construct the equations from the rightmost point
494  if ((*this->l_yend)[i]==0.0) {
495  ty[j]=yend2[i];
496  } else {
497  ty[j]=((*this->l_yend)[i]-yend2[i])/(*this->l_yend)[i];
498  }
499  j++;
500  }
501  }
502 
503  return 0;
504  }
505 
506 #endif
507 
508  };
509 
510 #ifndef DOXYGEN_NO_O2NS
511 }
512 #endif
513 
514 #endif
static const int right
Known on the right boundary.
Definition: ode_bv_solve.h:59
int set_mroot(mroot< mm_funct<> > &root)
Set the equation solver.
Definition: ode_bv_solve.h:237
double l_h
Storage for the stepsize.
Definition: ode_bv_solve.h:280
Solve an initial-value ODE problems given an adaptive ODE stepper.
Definition: ode_iv_solve.h:91
gsl_mroot_hybrids< mm_funct<> > def_mroot
The default equation solver.
Definition: ode_bv_solve.h:246
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
size_t mem_size
Size of recent allocation.
Definition: ode_bv_solve.h:294
int solve_grid_fun(size_t nv, const vec_t &tx, vec_t &ty)
The shooting function to be solved by the multidimensional solver.
Definition: ode_bv_solve.h:468
int solve_final_value(double x0, double x1, double h, size_t n, vec_t &ystart, vec_t &yend, vec_int_t &index, vec_t &yerr, vec_t &dydx_end, func_t &derivs)
Solve the boundary-value problem and store the solution.
Definition: ode_bv_solve.h:124
vec_int_t * l_index
The index defining the boundary conditions.
Definition: ode_bv_solve.h:259
int verbose
Set output level.
Definition: ode_bv_solve.h:67
static const int both
Known on both the left and right boundaries.
Definition: ode_bv_solve.h:63
ode_iv_solve< func_t, vec_t > def_ois
The default initial value solver.
Definition: ode_bv_solve.h:243
Multidimensional root-finding [abstract base].
Definition: mroot.h:66
Solve boundary-value ODE problems by shooting from one boundary to the other on a grid...
Definition: ode_bv_solve.h:359
static const int left
Known on the left boundary.
Definition: ode_bv_solve.h:61
vec_t * l_ystart
Storage for the starting vector.
Definition: ode_bv_solve.h:262
func_t * l_derivs
The functions to integrate.
Definition: ode_bv_solve.h:283
void allocate(size_t n)
Allocate internal storage.
Definition: ode_bv_solve.h:101
Solve boundary-value ODE problems by shooting from one boundary to the other.
Definition: ode_bv_solve.h:84
One-dimensional solver [abstract base].
Definition: root.h:47
vec_t * l_yend
Storage for the ending vector.
Definition: ode_bv_solve.h:265
mroot< mm_funct<> > * mrootp
The equation solver.
Definition: ode_bv_solve.h:256
int solve_final_value(double x0, double x1, double h, size_t n, vec_t &ystart, vec_t &yend, func_t &derivs)
Solve the initial-value problem to get the final value.
Definition: ode_iv_solve.h:181
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
int solve_grid(double x0, double x1, double h, size_t n, vec_t &ystart, vec_t &yend, vec_int_t &index, size_t nsol, vec_t &xsol, mat_t &ysol, mat_t &err_sol, mat_t &dydx_sol, func_t &derivs)
Desc.
Definition: ode_bv_solve.h:375
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
int solve_fun(size_t nv, const vec_t &tx, vec_t &ty)
The shooting function to be solved by the multidimensional solver.
Definition: ode_bv_solve.h:297
int set_iv(ode_iv_solve< func_t, vec_t > &ois)
Set initial value solver.
Definition: ode_bv_solve.h:231
size_t l_n
The number of functions.
Definition: ode_bv_solve.h:286
vec_t * l_yerr
Storage for the starting vector.
Definition: ode_bv_solve.h:268
static const int unk
Unknown on both the left and right boundaries.
Definition: ode_bv_solve.h:57
ode_iv_solve< func_t, vec_t > * oisp
The solver for the initial value problem.
Definition: ode_bv_solve.h:253
double l_x1
Storage for the ending abcissa.
Definition: ode_bv_solve.h:277
vec_t * l_dydx_end
Storage for the ending vector.
Definition: ode_bv_solve.h:271
double l_x0
Storage for the starting point.
Definition: ode_bv_solve.h:274
static const double x1[5]
Definition: inte_qng_gsl.h:48
Base class for boundary-value ODE solvers.
Definition: ode_bv_solve.h:44
std::function< int(double, size_t, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &)> ode_funct11
Ordinary differential equation function.
Definition: ode_funct.h:46
int solve_store(double x0, double x1, double h, size_t n, vec_t &ystart, vec_t &yend, vec_int_t &index, size_t &n_sol, vec_t &x_sol, mat_t &y_sol, mat_t &yerr_sol, mat_t &dydx_sol, func_t &derivs)
Solve the boundary-value problem and store the solution.
Definition: ode_bv_solve.h:196

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).