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

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