ode_bv_mshoot.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2008-2014, Julien Garaud and
5  Andrew W. Steiner
6 
7  This file is part of O2scl.
8 
9  O2scl is free software; you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation; either version 3 of the License, or
12  (at your option) any later version.
13 
14  O2scl is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  GNU General Public License for more details.
18 
19  You should have received a copy of the GNU General Public License
20  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
21 
22  -------------------------------------------------------------------
23 */
24 #ifndef O2SCL_ODE_BV_MSHOOT_H
25 #define O2SCL_ODE_BV_MSHOOT_H
26 
27 /** \file ode_bv_mshoot.h
28  \brief File defining \ref o2scl::ode_bv_mshoot
29 */
30 
31 #include <o2scl/ode_bv_solve.h>
32 
33 #ifndef DOXYGEN_NO_O2NS
34 namespace o2scl {
35 #endif
36 
37  /** \brief Solve boundary-value ODE problems by multishooting
38  with a generic nonlinear solver
39 
40  This class is experimental.
41 
42  Default template arguments
43  - \c func_t - \ref ode_funct11
44  - \c mat_t - \ref boost::numeric::ublas::matrix < double >
45  - \c vec_t - \ref boost::numeric::ublas::vector < double >
46  - \c vec_int_t - \ref boost::numeric::ublas::vector < int >
47 
48  \future Make a class which performs an iterative linear
49  solver which uses sparse matrices like ode_it_solve?
50  */
51  template<class func_t=ode_funct11,
54  class vec_int_t=boost::numeric::ublas::vector<int> >
55  class ode_bv_mshoot : public ode_bv_solve {
56 
57  public:
58 
59  ode_bv_mshoot() {
60  oisp=&def_ois;
62  mem_size=0;
63  }
64 
65  virtual ~ode_bv_mshoot() {
66  }
67 
68  /** \brief Solve the boundary-value problem and store the
69  solution
70  */
71  int solve_final_value(double h, size_t n, size_t n_bound,
72  vec_t &x_bound, mat_t &y_bound,
73  vec_int_t &index, func_t &derivs) {
74 
75  // Make copies of the inputs for later access
76  this->l_index=&index;
77  this->l_nbound=n_bound;
78  this->l_xbound=&x_bound;
79  this->l_ybound=&y_bound;
80  this->l_h=h;
81  this->l_derivs=&derivs;
82  this->l_n=n;
83 
84  int lhs_unks=0, rhs_conts=0;
85  this->l_lhs_unks=lhs_unks;
86 
87  // Count the number of variables we'll need to solve for
88  for (size_t i=0;i<n;i++) {
89  if (index[i]<2) lhs_unks++;
90  if (index[i]%2==1) rhs_conts++;
91  }
92 
93  // Make sure that the boundary conditions make sense
94  if (lhs_unks!=rhs_conts) {
95  O2SCL_ERR2("Incorrect boundary conditions in ",
96  "ode_bv_mshoot::solve()",gsl_einval);
97  }
98 
99  // The number of variables to solve for
100  int n_solve=lhs_unks+n*(n_bound-2);
101 
102  // Make space for the solution
103  ubvector tx(n_solve);
104 
105  // Copy initial guess from y_bound
106  size_t j=0;
107  for (size_t i=0;i<n;i++) {
108  if (index[i]<2) {
109  tx[j]=y_bound[0][i];
110  j++;
111  }
112  }
113  for(size_t k=1;k<n_bound-1;k++) {
114  for(size_t i=0;i<n;i++) {
115  tx[j]=y_bound[k][i];
116  j++;
117  }
118  }
119 
120  // Allocate memory as needed
121  if (n!=mem_size) {
122  sy.resize(n);
123  sy2.resize(n);
124  syerr.resize(n);
125  sdydx.resize(n);
126  sdydx2.resize(n);
127  mem_size=n;
128  }
129 
130  // The function object
131  mm_funct_mfptr<ode_bv_mshoot<func_t,mat_t,vec_t,vec_int_t> >
133 
134  // Solve
135  int ret=this->mrootp->msolve(n_solve,tx,mfm);
136  if (ret!=0) {
137  O2SCL_ERR("Solver failed in ode_bv_mshoot::solve().",ret);
138  }
139 
140  return ret;
141  }
142 
143  /** \brief Solve the boundary-value problem and store the
144  solution
145  */
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) {
151 
152  if (n_bound<2) {
153  O2SCL_ERR2("Not enough boundaries (must be at least two) in ",
154  "ode_bv_mshoot::solve_store().",gsl_einval);
155  }
156  if (n_sol<n_bound) {
157  O2SCL_ERR2("Not enough room to store boundaries in ",
158  "ode_bv_mshoot::solve_store().",gsl_einval);
159  }
160 
161  // Solve to fix x_bound and y_bound
162  solve_final_value(h,n,n_bound,x_bound,y_bound,index,derivs);
163 
164  // Find the indices which correspond to the boundaries
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;
170  }
171  // Double check that each interval has some space
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);
176  }
177  }
178 
179  // Now create the table
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);
183 
184  std::cout << "Old boundaries: " << inxs << std::endl;
185 
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]);
189 
190  std::cout << "New boundaries: " << n_sol_tmp << std::endl;
191 
192  // If it didn't use all the space, shift the indexes
193  // accordingly
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;
198  }
199  }
200 
201  std::cout << "New boundaries: " << inxs << std::endl;
202  }
203 
204  n_sol=inxs[n_bound-1];
205 
206  return 0;
207  }
208 
209  /** \brief Set initial value solver
210  */
212  oisp=&ois;
213  return 0;
214  }
215 
216  /** \brief Set the equation solver */
217  int set_mroot(mroot<mm_funct<> > &root) {
218  mrootp=&root;
219  return 0;
220  }
221 
222  /// The default initial value solver
224 
225  /// The default equation solver
226  gsl_mroot_hybrids<mm_funct<> > def_mroot;
227 
228 #ifndef DOXYGEN_INTERNAL
229 
230  protected:
231 
232  /// The solver for the initial value problem
234 
235  /// The equation solver
237 
238  /// The index defining the boundary conditions
239  vec_int_t *l_index;
240 
241  /// Storage for the starting vector
242  vec_t *l_xbound;
243 
244  /// Storage for the ending vector
245  mat_t *l_ybound;
246 
247  /// Storage for the stepsize
248  double l_h;
249 
250  /// The functions to integrate
251  func_t *l_derivs;
252 
253  /// The number of functions
254  size_t l_n;
255 
256  /// The number of boundaries
257  size_t l_nbound;
258 
259  /// The number of unknowns on the LHS
260  size_t l_lhs_unks;
261 
262  /// \name Temporary storage for \ref solve_fun()
263  //@{
264  vec_t sy, sy2, syerr, sdydx, sdydx2;
265  //@}
266 
267  /// Size of recent allocation
268  size_t mem_size;
269 
270  /// The shooting function to be solved by the multidimensional solver
271  int solve_fun(size_t nv, const vec_t &tx, vec_t &ty) {
272 
273  // Counters to index through function parameters tx and ty
274  size_t j_x=0, j_y=0;
275 
276  // Shorthand for the boundary specifications
277  vec_t &xb=*this->l_xbound;
278  mat_t &yb=*this->l_ybound;
279 
280  // Set up the boundaries from the values in tx
281  for(size_t k=0;k<l_nbound-1;k++) {
282  if (k==0) {
283  for(size_t i=0;i<this->l_n;i++) {
284  if ((*this->l_index)[i]<2) {
285  yb[k][i]=tx[j_x];
286  j_x++;
287  }
288  }
289  } else {
290  for(size_t i=0;i<this->l_n;i++) {
291  yb[k][i]=tx[j_x];
292  j_x++;
293  }
294  }
295  }
296 
297  // Integrate between all of the boundaries
298  for(size_t k=0;k<l_nbound-1;k++) {
299 
300  double x0=xb[k];
301  double x1=xb[k+1];
302 
303  // Setup the start vector sy. This extra copying from l_ybound
304  // to sy might be avoided by using a mat_row_t object to send
305  // l_ybound directly to the solver?
306 
307  for(size_t i=0;i<this->l_n;i++) {
308  sy[i]=yb[k][i];
309  }
310 
311  // Shoot across. We specify memory for the derivatives and
312  // errors to prevent ode_iv_solve from having to continuously
313  // allocate and reallocate it.
314  oisp->solve_final_value(x0,x1,this->l_h,this->l_n,sy,
315  sy2,*this->l_derivs);
316 
317  if (k!=l_nbound-2) {
318  // Construct equations for the internal RHS boundaries
319  for(size_t i=0;i<this->l_n;i++) {
320  ty[j_y]=sy2[i]-yb[k+1][i];
321  j_y++;
322  }
323  } else {
324  // Construct the equations for the rightmost boundary
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];
328  if (yright==0.0) {
329  ty[j_y]=sy2[i];
330  } else {
331  ty[j_y]=(sy2[i]-yright)/yright;
332  }
333  j_y++;
334  }
335  }
336  }
337 
338  // End of loop to integrate between all boundaries
339  }
340 
341  return 0;
342  }
343 
344 #endif
345 
346  };
347 
348 #ifndef DOXYGEN_NO_O2NS
349 }
350 #endif
351 
352 #endif
Solve an initial-value ODE problems given an adaptive ODE stepper.
Definition: ode_iv_solve.h:91
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=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
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.
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].
Definition: mroot.h:66
Solve boundary-value ODE problems by multishooting with a generic nonlinear solver.
Definition: ode_bv_mshoot.h:55
size_t l_nbound
The number of boundaries.
func_t * l_derivs
The functions to integrate.
One-dimensional solver [abstract base].
Definition: root.h:47
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 &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_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.
Definition: ode_bv_mshoot.h:71
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.
Definition: err_hnd.h:273
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]
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

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