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

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