ode_iv_solve.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, 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_IV_SOLVE_H
24 #define O2SCL_ODE_IV_SOLVE_H
25 
26 /** \file ode_iv_solve.h
27  \brief File defining \ref o2scl::ode_iv_solve
28 */
29 
30 #include <string>
31 
32 #include <boost/numeric/ublas/vector.hpp>
33 #include <boost/numeric/ublas/matrix.hpp>
34 #include <boost/numeric/ublas/matrix_proxy.hpp>
35 
36 #include <o2scl/astep.h>
37 #include <o2scl/astep_gsl.h>
38 
39 #ifndef DOXYGEN_NO_O2NS
40 namespace o2scl {
41 #endif
42 
43  /** \brief Solve an initial-value ODE problems given an adaptive ODE
44  stepper
45 
46  This class gives several functions which solve an initial
47  value ODE problem. The functions \ref solve_final_value()
48  gives only the final value of the functions at the end
49  of the ODE integration and is relatively fast.
50 
51  The function solve_store() stores the solution of the ODE over
52  the full range into a set of vectors and matrices which are
53  allocated and specified by the user. This function is designed
54  to give exactly the same results (though this cannot be
55  guaranteed) as solve_final_value() and additionally records some
56  or all of the results from the adaptive steps which were taken.
57 
58  All of these functions automatically evaluate the derivatives
59  from the specified function at the initial point and
60  user-specified initial derivatives are ignored. The total
61  number of steps taken is limited by \ref ntrial and \ref
62  nsteps stores the number of steps taken by the most recent
63  solution. The variable \ref nsteps_out is the maximum number
64  of points in the interval for which verbose output will be
65  given when \ref o2scl::ode_iv_solve::verbose is greater than zero.
66 
67  There is an example for the usage of this class in
68  <tt>examples/ex_ode.cpp</tt> documented in the \ref ex_ode_sect
69  section.
70 
71  <b>Convergence error handling</b>
72 
73  There are two different convergence errors which can
74  be controlled separately in this class.
75  - The adaptive stepper may require too many steps. If this
76  happens, then the solver immediately stops. The solver
77  calls the error handler if \ref err_nonconv is true, and
78  otherwise it returns a non-zero value.
79  - The adaptive stepper may fail. If \ref exit_on_fail
80  is true, then the error handler is called. Otherwise,
81  the solver proceeds to continue computing the whole solution.
82  So long as the number of adaptive steps required is less
83  than \ref ntrial, then the full solution is computed and
84  a non-zero value is returned to indicate the accuracy of
85  the solution may be impacted. If the number of adaptive
86  steps required after a failure of the adaptive stepper
87  is larger than \ref ntrial, then the behavior of the
88  solver is controlled by \ref err_nonconv as described
89  above.
90 
91  Documentation links for default template arguments
92  - \c func_t - \ref ode_funct
93  - \c vec_t - \ref boost::numeric::ublas::vector < double >
94 
95  The default adaptive stepper is an object of type \ref astep_gsl.
96 
97  \future The form of solve_final_value() is very similar to that
98  of astep_base::astep_full(), but not quite the same. Maybe
99  these functions should be consistent with each other?
100  */
101  template<class func_t=ode_funct,
103  class ode_iv_solve {
104 
105  public:
106 
108 
109 #ifndef DOXYGEN_INTERNAL
110 
111  protected:
112 
113  /// \name Vectors for temporary storage
114  //@{
115  vec_t vtemp, vtemp2, vtemp3, vtemp4;
116  //@}
117 
118  /// The size of the temporary vectors
119  size_t mem_size;
120 
121  /// The adaptive stepper
123 
124  /// Print out iteration information
125  virtual int print_iter(double x, size_t nv, vec_t &y) {
126  std::cout << type() << " x: " << x << " y: ";
127  for(size_t i=0;i<nv;i++) std::cout << y[i] << " ";
128  std::cout << std::endl;
129  if (verbose>1) {
130  char ch;
131  std::cin >> ch;
132  }
133  return 0;
134  }
135 
136  /// Free allocated memory
137  void free() {
138  if (mem_size>0) {
139  vtemp.clear();
140  vtemp2.clear();
141  vtemp3.clear();
142  vtemp4.clear();
143  }
144  }
145 
146  /** \brief Allocate space for temporary vectors
147  */
148  void allocate(size_t n) {
149  if (n!=mem_size) {
150  free();
151  vtemp.resize(n);
152  vtemp2.resize(n);
153  vtemp3.resize(n);
154  vtemp4.resize(n);
155  mem_size=n;
156  }
157  }
158 
159 #endif
160 
161  public:
162 
163  ode_iv_solve() {
164  verbose=0;
165  ntrial=1000;
166  nsteps_out=10;
167  astp=&gsl_astp;
168  exit_on_fail=true;
169  mem_size=0;
170  err_nonconv=true;
171  }
172 
173  virtual ~ode_iv_solve() {
174  free();
175  }
176 
177  /** \brief If true, call the error handler if the solution does
178  not converge (default true)
179  */
181 
182  /// \name Main solver functions
183  //@{
184  /** \brief Solve the initial-value problem to get the final value
185 
186  Given the \c n initial values of the functions in \c ystart,
187  this function integrates the ODEs specified in \c derivs over
188  the interval from \c x0 to \c x1 with an initial stepsize of \c
189  h. The final values of the function are given in \c yend and the
190  initial values of \c yend are ignored.
191 
192  If \ref verbose is greater than zero, The solution at less than
193  or approximately equal to \ref nsteps_out points will be written
194  to \c std::cout. If \ref verbose is greater than one, a
195  character will be required after each selected point.
196  */
197  int solve_final_value(double x0, double x1, double h, size_t n,
198  vec_t &ystart, vec_t &yend, func_t &derivs) {
199 
200  allocate(n);
201  return solve_final_value(x0,x1,h,n,ystart,yend,vtemp2,
202  vtemp3,derivs);
203  }
204 
205  /** \brief Solve the initial-value problem to get the final value
206  with errors
207 
208  Given the \c n initial values of the functions in \c ystart,
209  this function integrates the ODEs specified in \c derivs over
210  the interval from \c x0 to \c x1 with an initial stepsize of \c
211  h. The final values of the function are given in \c yend and the
212  associated errors are given in \c yerr. The initial values of \c
213  yend and \c yerr are ignored.
214 
215  If \ref verbose is greater than zero, The solution at less
216  than or approximately equal to \ref nsteps_out points will be
217  written to \c std::cout. If \ref verbose is greater than one,
218  a character will be required after each selected point.
219  */
220  int solve_final_value(double x0, double x1, double h, size_t n,
221  vec_t &ystart, vec_t &yend, vec_t &yerr,
222  func_t &derivs) {
223 
224  allocate(n);
225  return solve_final_value(x0,x1,h,n,ystart,yend,yerr,
226  vtemp2,derivs);
227  }
228 
229  /** \brief Solve the initial-value problem to get the final value,
230  derivatives, and errors
231 
232  Given the \c n initial values of the functions in \c ystart,
233  this function integrates the ODEs specified in \c derivs over
234  the interval from \c x0 to \c x1 with an initial stepsize of \c
235  h. The final values of the function are given in \c yend, the
236  derivatives in \c dydx_end, and the associated errors are given
237  in \c yerr. The initial values of \c yend and \c yerr are
238  ignored.
239 
240  This function is designed to be relatively fast,
241  avoiding extra copying of vectors back and forth.
242 
243  If \ref verbose is greater than zero, The solution at less
244  than or approximately equal to \ref nsteps_out points will be
245  written to \c std::cout. If \ref verbose is greater than one,
246  a character will be required after each selected point.
247 
248  This function computes \c dydx_start automatically and the
249  values given by the user are ignored.
250 
251  The solution fails if more than \ref ntrial steps are required.
252  This function will also fail if <tt>x1>x0</tt> and <tt>h<0</tt>
253  or if <tt>x1-x0</tt> and <tt>h>0</tt> do not have the same sign.
254  */
255  int solve_final_value(double x0, double x1, double h, size_t n,
256  vec_t &ystart, vec_t &yend, vec_t &yerr,
257  vec_t &dydx_end, func_t &derivs) {
258 
259  if ((x1>x0 && h<=0.0) || (x0>x1 && h>=0.0)) {
260  std::string str="Interval direction (x1-x0="+o2scl::dtos(x1-x0)+
261  ") does not match step direction (h="+o2scl::dtos(h)+
262  " in ode_iv_solve::solve_final_value().";
263  O2SCL_ERR(str.c_str(),exc_einval);
264  }
265  if (x0==x1) {
266  O2SCL_ERR2("Starting and final endpoints identical in ",
267  "ode_iv_solve::solve_final_value().",exc_einval);
268  }
269 
270  // Allocate for temporary vectors
271  allocate(n);
272 
273  // The variable 'x' is the current independent variable, xnext is
274  // the next value of x, xverb is the next value of x for which
275  // verbose output should be provided, and dxverb is the stepsize
276  // for xverb.
277  double x=x0, xverb=0.0, dxverb=0.0, xnext;
278  int ret=0, first_ret=0;
279 
280  nsteps=0;
281 
282  // Create a reference to make the code easier to read
283  vec_t &dydx_start=vtemp;
284 
285  // Compute stepsize for verbose output
286  if (verbose>0) {
287  print_iter(x0,n,ystart);
288  if (verbose>1) {
289  char ch;
290  std::cin >> ch;
291  }
292  // Ensure that 'dxverb' is positive
293  if (x1>x0) {
294  dxverb=(x1-x0)/((double)nsteps_out);
295  xverb=x0+dxverb;
296  } else {
297  dxverb=(x0-x1)/((double)nsteps_out);
298  xverb=x0-dxverb;
299  }
300  }
301 
302  // Use yend as workspace
303  for(size_t i=0;i<n;i++) yend[i]=ystart[i];
304 
305  // Compute initial derivative
306  derivs(x,n,yend,dydx_start);
307 
308  bool done=false;
309  while (done==false) {
310 
311  // Take a step of the first type
312  ret=astp->astep_full(x,x1,xnext,h,n,ystart,dydx_start,
313  yend,yerr,dydx_end,derivs);
314 
315  if (ret!=0) {
316  if (exit_on_fail) {
317  O2SCL_ERR2("Adaptive stepper failed in ",
318  "ode_iv_solve::solve_final_value()",ret);
319  } else if (first_ret!=0) {
320  first_ret=ret;
321  }
322  }
323 
324  // Print out verbose info
325  if (verbose>0 && xnext>xverb) {
326  print_iter(xnext,n,yend);
327  if (verbose>1) {
328  char ch;
329  std::cin >> ch;
330  }
331  xverb+=dxverb;
332  }
333 
334  // Check number of iterations
335  nsteps++;
336  if (nsteps>ntrial) {
337  std::string str="Too many steps required (nsteps="+
338  szttos(nsteps)+", ntrial="+szttos(ntrial)+
339  ", x="+o2scl::dtos(x)+") in ode_iv_solve::solve_final_value().";
341  }
342 
343  if (ret!=0) {
344  done=true;
345  } else {
346  if (x1>x0) {
347  if (xnext>=x1) done=true;
348  } else {
349  if (xnext<=x1) done=true;
350  }
351  }
352 
353  if (done==false) {
354 
355  // Take a step of the second type
356  ret=astp->astep_full(xnext,x1,x,h,n,yend,dydx_end,ystart,yerr,
357  dydx_start,derivs);
358 
359  if (ret!=0) {
360  if (exit_on_fail) {
361  O2SCL_ERR2("Adaptive stepper failed in ",
362  "ode_iv_solve::solve_final_value()",ret);
363  } else if (first_ret!=0) {
364  first_ret=ret;
365  }
366  }
367 
368  // Print out verbose info
369  if (verbose>0 && x>xverb) {
370  print_iter(x,n,ystart);
371  if (verbose>1) {
372  char ch;
373  std::cin >> ch;
374  }
375  xverb+=dxverb;
376  }
377 
378  // Check number of iterations
379  nsteps++;
380  if (nsteps>ntrial) {
381  std::string str="Too many steps required (ntrial="+itos(ntrial)+
382  ", x="+o2scl::dtos(x)+") in ode_iv_solve::solve_final_value().";
383  O2SCL_ERR(str.c_str(),exc_emaxiter);
384  }
385 
386  if (ret!=0) {
387  done=true;
388  } else {
389  if (x1>x0) {
390  if (x>=x1) {
391  done=true;
392  for(size_t j=0;j<n;j++) {
393  yend[j]=ystart[j];
394  dydx_end[j]=dydx_start[j];
395  }
396  }
397  } else {
398  if (x<=x1) {
399  done=true;
400  for(size_t j=0;j<n;j++) {
401  yend[j]=ystart[j];
402  dydx_end[j]=dydx_start[j];
403  }
404  }
405  }
406  }
407 
408  }
409 
410  // End of while loop
411  }
412 
413  // Print out final verbose info
414  if (verbose>0) {
415  print_iter(x1,n,yend);
416  if (verbose>1) {
417  char ch;
418  std::cin >> ch;
419  }
420  }
421 
422  return first_ret;
423  }
424 
425  /** \brief Solve the initial-value problem and store the
426  associated output
427 
428  Initially, \c x_sol should be a vector of size \c n_sol, and \c
429  y_sol, \c dydx_sol, and \c yerr_sol should be matrices with size
430  \c [n_sol][n]. On exit, \c n_sol will will be number of points
431  store, less than or equal to the original value of \c
432  n_sol. This function avoids performing extra calls to the
433  adaptive stepper, and the solution will be approximately evenly
434  spaced.
435 
436  This function is also designed to give the exactly the same
437  results as solve_final_value(). This cannot be strictly
438  guaranteed, but will likely hold in most applications.
439 
440  This template function works with any matrix class \c mat_t
441  which can be accessed using successive applications of
442  operator[] and which has an associated class \c mat_row_t
443  which returns a row of a matrix of type \c mat_t as
444  an object with type \c vec_t.
445 
446  If \ref verbose is greater than zero, The solution at each
447  internal point will be written to \c std::cout. If \ref verbose
448  is greater than one, a character will be required after each
449  point.
450 
451  \todo Document \c istart
452  */
453  template<class mat_t>
454  int solve_store(double x0, double x1, double h, size_t n, size_t &n_sol,
455  vec_t &x_sol, mat_t &y_sol, mat_t &yerr_sol,
456  mat_t &dydx_sol, func_t &derivs, size_t istart=0) {
457 
458  int ret=0;
459  int first_ret=0;
460  size_t nmax=n_sol-1;
461  nsteps=0;
462 
463  // Stepsize for next verbose output. Use nsteps_out-1 instead of
464  // nsteps_out since the first point is always output below.
465  double dx_verb=(x1-x0)/((double)(nsteps_out-1));
466  // Stepsize for next point for storage
467  double dx_tab=(x1-x0)/((double)(n_sol-istart-1));
468 
469  double x_verb=x0+dx_verb;
470  double x_tab=x0+dx_tab;
471 
472  // Allocate space for errors and derivatives and extra storage
473  allocate(n);
474 
475  // Create some references just to make the code easier
476  // to read
477  vec_t &ystart=vtemp4;
478  vec_t &dydx_start=vtemp2;
479 
480  // Copy first point to ystart
481  for(size_t j=0;j<n;j++) ystart[j]=y_sol(istart,j);
482 
483  // Output first point
484  if (verbose>0) {
485  print_iter(x0,n,ystart);
486  if (verbose>1) {
487  char ch;
488  std::cin >> ch;
489  }
490  }
491 
492  // Initial derivative evaulation
493  derivs(x0,n,ystart,dydx_start);
494 
495  // Add first derivatives to storage, and set the errors to zero
496  x_sol[istart]=x0;
497  for(size_t j=0;j<n;j++) {
498  dydx_sol(istart,j)=dydx_start[j];
499  yerr_sol(istart,j)=0.0;
500  }
501 
502  // Copy first point to storage again for first step
503  size_t icurr=istart+1;
504  x_sol[icurr]=x0;
505  for(size_t j=0;j<n;j++) y_sol(icurr,j)=ystart[j];
506 
507  double xnext;
508 
509  bool done=false;
510  while (done==false) {
511 
512  // Create some references just to make the code easier
513  // to read
514  vec_t &yerr=vtemp;
515  vec_t &dydx_out=vtemp3;
516 
517  // Use ystart as temporary storage for the end of the current
518  // adaptive step
519  vec_t yrow(n);
520  for(size_t i=0;i<n;i++) yrow[i]=y_sol(icurr,i);
521  ret=astp->astep_full(x0,x1,xnext,h,n,yrow,dydx_start,ystart,yerr,
522  dydx_out,derivs);
523 
524  if (ret!=0) {
525  if (exit_on_fail) {
526  n_sol=icurr+1;
527  // call error handler
528  O2SCL_ERR2("Adaptive stepper returned non-zero in ",
529  "ode_iv_solve::solve_store().",exc_efailed);
530  } else if (first_ret==0) {
531  first_ret=ret;
532  }
533  }
534 
535  // Update step counter and abort if necessary
536  nsteps++;
537  if (nsteps>ntrial) {
538  std::string str="Too many steps required (ntrial="+itos(ntrial)+
539  ", x="+o2scl::dtos(x0)+") in ode_iv_solve::solve_store().";
540  O2SCL_ERR(str.c_str(),exc_emaxiter);
541  }
542 
543  // If we've made enough progress, do verbose output
544  if (xnext>=x_verb && verbose>0) {
545  print_iter(xnext,n,ystart);
546  if (verbose>1) {
547  char ch;
548  std::cin >> ch;
549  }
550  x_verb+=dx_verb;
551  }
552 
553  // If we're at the end
554  if (xnext>=x1) {
555 
556  // Exit the loop
557  done=true;
558 
559  // Store the final entry
560  x_sol[icurr]=xnext;
561  for(size_t j=0;j<n;j++) {
562  y_sol(icurr,j)=ystart[j];
563  dydx_sol(icurr,j)=dydx_out[j];
564  yerr_sol(icurr,j)=yerr[j];
565  }
566 
567  // Update the solution size
568  n_sol=icurr+1;
569 
570  } else {
571 
572  if (xnext>=x_tab && icurr<nmax) {
573 
574  // If we've made enough progress, store an entry in the table
575  x_sol[icurr]=xnext;
576  for(size_t j=0;j<n;j++) {
577  y_sol(icurr,j)=ystart[j];
578  dydx_sol(icurr,j)=dydx_out[j];
579  yerr_sol(icurr,j)=yerr[j];
580 
581  // Also prepare for the next step
582  y_sol(icurr+1,j)=ystart[j];
583  dydx_start[j]=dydx_out[j];
584  }
585  x_tab+=dx_tab;
586 
587  // Update x0 and the current and next indicies
588  x0=xnext;
589  icurr++;
590 
591  } else {
592 
593  // Otherwise, just prepare for the next step
594  // without updating icurr
595  x0=xnext;
596  for(size_t j=0;j<n;j++) {
597  dydx_start[j]=dydx_out[j];
598  y_sol(icurr,j)=ystart[j];
599  }
600 
601  }
602 
603  }
604 
605  // End of loop 'while (done==false)'
606  }
607 
608  return first_ret;
609  }
610  //@}
611 
612  /// Set output level
613  int verbose;
614 
615  /** \brief Number of output points for verbose output (default 10)
616 
617  This is used in functions solve_store() and solve_final_value()
618  to control how often steps are output when verbose is greater
619  than zero.
620  */
621  size_t nsteps_out;
622 
623  /// Maximum number of applications of the adaptive stepper (default 1000)
624  size_t ntrial;
625 
626  /// Number of adaptive ste!ps employed
627  size_t nsteps;
628 
629  /// \name The adaptive stepper
630  //@{
631  /// Set the adaptive stepper to use
633  astp=&as;
634  return 0;
635  }
636 
637  /** \brief If true, stop the solution if the adaptive stepper fails
638  (default true)
639  */
641 
642  /// The default adaptive stepper
644  //@}
645 
646  /// Return the type, \c "ode_iv_solve".
647  virtual const char *type() { return "ode_iv_solve"; }
648 
649  };
650 
651  typedef matrix_row_gen_ctor<boost::numeric::ublas::matrix<double> >
652  solve_grid_mat_row;
653 
654  typedef std::function<int(double,size_t,const solve_grid_mat_row &,
655  solve_grid_mat_row &)>
656  ode_funct_solve_grid;
657 
658  /** \brief Solve an initial-value ODE problems on a grid
659  given an adaptive ODE stepper
660 
661  This class works as similar to ode_iv_solve::solve_store()
662  except that the solution is stored on a grid of points in the
663  independent variable specified by the user, at the cost of
664  taking extra steps to ensure that function values, derivatives,
665  and errors are computed at each grid point.
666 
667  There is an example for the usage of this class in
668  <tt>examples/ex_ode.cpp</tt> documented in the \ref ex_ode_sect
669  section.
670  */
671  template<class func_t=ode_funct_solve_grid,
672  class mat_row_t=solve_grid_mat_row>
674 
675 #ifndef DOXYGEN_INTERNAL
676 
677  protected:
678 
679  /// The adaptive stepper
681 
682  /// Print out iteration information
683  virtual int print_iter(double x, size_t nv, mat_row_t &y) {
684  std::cout << type() << " x: " << x << " y: ";
685  for(size_t i=0;i<nv;i++) std::cout << y[i] << " ";
686  std::cout << std::endl;
687  if (verbose>1) {
688  char ch;
689  std::cin >> ch;
690  }
691  return 0;
692  }
693 
694 #endif
695 
696  public:
697 
699  verbose=0;
700  ntrial=1000;
701  astp=&gsl_astp;
702  exit_on_fail=true;
703  err_nonconv=true;
704  }
705 
706  virtual ~ode_iv_solve_grid() {
707  }
708 
709  /** \brief If true, call the error handler if the solution does
710  not converge (default true)
711  */
713 
714  /// \name Main solver function
715  //@{
716  /** \brief Solve the initial-value problem from \c x0 to \c x1 over
717  a grid storing derivatives and errors
718 
719  Initially, \c xsol should be an array of size \c nsol, and \c
720  ysol should be a \c ubmatrix of size \c [nsol][n]. This function
721  never takes a step larger than the grid size.
722 
723  If \ref verbose is greater than zero, The solution at each grid
724  point will be written to \c std::cout. If \ref verbose is
725  greater than one, a character will be required after each point.
726 
727  \future Consider making a version of grid which gives the
728  same answers as solve_final_value(). After each proposed step,
729  it would go back and fill in the grid points if the proposed
730  step was past the next grid point.
731  */
732  template<class vec_t, class mat_t>
733  int solve_grid(double h, size_t n, size_t nsol, vec_t &xsol,
734  mat_t &ysol, mat_t &err_sol, mat_t &dydx_sol,
735  func_t &derivs) {
736 
737  double x0=xsol[0];
738  double x1=xsol[nsol-1];
739 
740  double x=x0, xnext;
741  int ret=0, first_ret=0;
742  nsteps=0;
743 
744  // Copy the initial point to the first row
745  xsol[0]=x0;
746 
747  mat_row_t y_start(n);
748  mat_row_t dydx_start(n);
749 
750  // Copy ysol[0] to ystart
751  for(size_t j=0;j<n;j++) {
752  y_start[j]=ysol(0,j);
753  }
754 
755  // Verbose output for the first row
756  if (verbose>0) print_iter(xsol[0],n,y_start);
757 
758  // Evaluate the first derivative
759  derivs(x0,n,y_start,dydx_start);
760 
761  // Set the derivatives and the uncertainties for the first row
762  for(size_t j=0;j<n;j++) {
763  dydx_sol(0,j)=dydx_start[j];
764  err_sol(0,j)=0.0;
765  }
766 
767  for(size_t i=1;i<nsol && ret==0;i++) {
768 
769  mat_row_t y_row(ysol,i);
770  mat_row_t dydx_row(dydx_sol,i);
771  mat_row_t yerr_row(err_sol,i);
772 
773  // Step until we reach the next grid point
774  bool done=false;
775  while(done==false && ret==0) {
776 
777  ret=astp->astep_full(x,xsol[i],xnext,h,n,y_start,dydx_start,
778  y_row,yerr_row,dydx_row,derivs);
779 
780  nsteps++;
781  if (ret!=0) {
782  if (exit_on_fail) {
783  O2SCL_ERR2("Adaptive stepper failed in ",
784  "ode_iv_solve_grid::solve_grid()",ret);
785  } else if (first_ret!=0) {
786  first_ret=ret;
787  }
788  }
789 
790  if (nsteps>ntrial) {
791  std::string str="Too many steps required (ntrial="+itos(ntrial)+
792  ", x="+o2scl::dtos(x)+") in ode_iv_solve_grid::solve_grid().";
793  O2SCL_ERR(str.c_str(),exc_emaxiter);
794  }
795 
796  // Adjust independent variable for next step
797  x=xnext;
798  for(size_t j=0;j<n;j++) {
799  y_start[j]=y_row[j];
800  dydx_start[j]=dydx_row[j];
801  }
802 
803  // Decide if we have reached the grid point
804  if (x1>x0) {
805  if (x>=xsol[i]) done=true;
806  } else {
807  if (x<=xsol[i]) done=true;
808  }
809 
810  }
811 
812  if (verbose>0) print_iter(xsol[i],n,y_row);
813 
814  }
815 
816  return first_ret;
817  }
818  //@}
819 
820  /// Set output level
821  int verbose;
822 
823  /// Maximum number of applications of the adaptive stepper (default 1000)
824  size_t ntrial;
825 
826  /// Number of adaptive steps employed
827  size_t nsteps;
828 
829  /// \name The adaptive stepper
830  //@{
831  /// Set the adaptive stepper to use
833  astp=&as;
834  return 0;
835  }
836 
837  /** \brief If true, stop the solution if the adaptive stepper fails
838  (default true)
839  */
841 
842  /// The default adaptive stepper
844  //@}
845 
846  /// Return the type, \c "ode_iv_solve".
847  virtual const char *type() { return "ode_iv_solve_grid"; }
848 
849  };
850 
851 #ifndef DOXYGEN_NO_O2NS
852 }
853 #endif
854 
855 #endif
o2scl::ode_iv_solve_grid::set_astep
int set_astep(astep_base< mat_row_t, mat_row_t, mat_row_t, func_t > &as)
Set the adaptive stepper to use.
Definition: ode_iv_solve.h:832
o2scl::ode_iv_solve
Solve an initial-value ODE problems given an adaptive ODE stepper.
Definition: ode_iv_solve.h:103
o2scl::ode_iv_solve::solve_final_value
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:197
o2scl::ode_iv_solve_grid::type
virtual const char * type()
Return the type, "ode_iv_solve".
Definition: ode_iv_solve.h:847
o2scl::ode_iv_solve::solve_store
int solve_store(double x0, double x1, double h, size_t n, size_t &n_sol, vec_t &x_sol, mat_t &y_sol, mat_t &yerr_sol, mat_t &dydx_sol, func_t &derivs, size_t istart=0)
Solve the initial-value problem and store the associated output.
Definition: ode_iv_solve.h:454
o2scl::dtos
std::string dtos(const fp_t &x, int prec=6, bool auto_prec=false)
Convert a double to a string.
Definition: string_conv.h:77
o2scl::ode_iv_solve::astp
astep_base< vec_t, vec_t, vec_t, func_t > * astp
The adaptive stepper.
Definition: ode_iv_solve.h:122
boost::numeric::ublas::vector< double >
o2scl::exc_efailed
@ exc_efailed
generic failure
Definition: err_hnd.h:61
o2scl::ode_iv_solve_grid::solve_grid
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.
Definition: ode_iv_solve.h:733
o2scl::ode_iv_solve_grid::ntrial
size_t ntrial
Maximum number of applications of the adaptive stepper (default 1000)
Definition: ode_iv_solve.h:824
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
o2scl::ode_iv_solve_grid::nsteps
size_t nsteps
Number of adaptive steps employed.
Definition: ode_iv_solve.h:827
o2scl::ode_iv_solve::gsl_astp
astep_gsl< vec_t, vec_t, vec_t, func_t > gsl_astp
The default adaptive stepper.
Definition: ode_iv_solve.h:643
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_iv_solve::set_astep
int set_astep(astep_base< vec_t, vec_t, vec_t, func_t > &as)
Set the adaptive stepper to use.
Definition: ode_iv_solve.h:632
o2scl::exc_emaxiter
@ exc_emaxiter
exceeded max number of iterations
Definition: err_hnd.h:73
o2scl::ode_iv_solve::solve_final_value
int solve_final_value(double x0, double x1, double h, size_t n, vec_t &ystart, vec_t &yend, vec_t &yerr, vec_t &dydx_end, func_t &derivs)
Solve the initial-value problem to get the final value, derivatives, and errors.
Definition: ode_iv_solve.h:255
o2scl::ode_iv_solve::solve_final_value
int solve_final_value(double x0, double x1, double h, size_t n, vec_t &ystart, vec_t &yend, vec_t &yerr, func_t &derivs)
Solve the initial-value problem to get the final value with errors.
Definition: ode_iv_solve.h:220
o2scl::ode_iv_solve::type
virtual const char * type()
Return the type, "ode_iv_solve".
Definition: ode_iv_solve.h:647
o2scl::ode_iv_solve_grid::print_iter
virtual int print_iter(double x, size_t nv, mat_row_t &y)
Print out iteration information.
Definition: ode_iv_solve.h:683
o2scl::ode_iv_solve::exit_on_fail
bool exit_on_fail
If true, stop the solution if the adaptive stepper fails (default true)
Definition: ode_iv_solve.h:640
o2scl::ode_iv_solve_grid::verbose
int verbose
Set output level.
Definition: ode_iv_solve.h:821
o2scl_inte_qng_coeffs::x1
static const double x1[5]
Definition: inte_qng_gsl.h:48
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
O2SCL_CONV_RET
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
Definition: err_hnd.h:297
o2scl::ode_iv_solve::mem_size
size_t mem_size
The size of the temporary vectors.
Definition: ode_iv_solve.h:119
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_iv_solve::ntrial
size_t ntrial
Maximum number of applications of the adaptive stepper (default 1000)
Definition: ode_iv_solve.h:624
o2scl::astep_gsl
Adaptive ODE stepper (GSL)
Definition: astep_gsl.h:268
o2scl::ode_iv_solve_grid::err_nonconv
bool err_nonconv
If true, call the error handler if the solution does not converge (default true)
Definition: ode_iv_solve.h:712
o2scl::ode_iv_solve::nsteps_out
size_t nsteps_out
Number of output points for verbose output (default 10)
Definition: ode_iv_solve.h:621
o2scl::ode_iv_solve::err_nonconv
bool err_nonconv
If true, call the error handler if the solution does not converge (default true)
Definition: ode_iv_solve.h:180
o2scl::ode_iv_solve::verbose
int verbose
Set output level.
Definition: ode_iv_solve.h:613
o2scl::itos
std::string itos(int x)
Convert an integer to a string.
o2scl::ode_iv_solve::print_iter
virtual int print_iter(double x, size_t nv, vec_t &y)
Print out iteration information.
Definition: ode_iv_solve.h:125
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::astep_base
Adaptive stepper [abstract base].
Definition: astep.h:49
o2scl::ode_iv_solve_grid::exit_on_fail
bool exit_on_fail
If true, stop the solution if the adaptive stepper fails (default true)
Definition: ode_iv_solve.h:840
o2scl::ode_iv_solve_grid::astp
astep_base< mat_row_t, mat_row_t, mat_row_t, func_t > * astp
The adaptive stepper.
Definition: ode_iv_solve.h:680
o2scl::ode_iv_solve::nsteps
size_t nsteps
Number of adaptive ste!ps employed.
Definition: ode_iv_solve.h:627
o2scl::ode_iv_solve::free
void free()
Free allocated memory.
Definition: ode_iv_solve.h:137
o2scl::ode_iv_solve::allocate
void allocate(size_t n)
Allocate space for temporary vectors.
Definition: ode_iv_solve.h:148
o2scl::ode_iv_solve_grid
Solve an initial-value ODE problems on a grid given an adaptive ODE stepper.
Definition: ode_iv_solve.h:673
o2scl::szttos
std::string szttos(size_t x)
Convert a size_t to a string.
o2scl::ode_iv_solve_grid::gsl_astp
astep_gsl< mat_row_t, mat_row_t, mat_row_t, func_t > gsl_astp
The default adaptive stepper.
Definition: ode_iv_solve.h:843

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